This online ressource is there to provide additional interactive visualization tool for the Varroa mites sampling and sequencing data obtained from the associated paper “XXXX”. A Snakemake pipeline created and used for genomics mapping, variant calling and analysis is provided in this host GitHub repository. Most genomics and demographic analysis using fastsimcoal2 were computed on the Okinawa Institute of Science and Technology OIST cluster, Sango.
Click on “code” button on the right to make appear the libraries used and R code for each section.
library("tidyverse")
library("knitr") # for the markdown
library("kableExtra") # for creating a scrolling table
library("ggplot2") # for plotting
library("maps")
library("leaflet")
library("htmltools")
library("rgdal")
library("gridExtra")
library("gghighlight")
library("adegenet")
library("vcfR")
library("ggmap")
library("poppr")
library("RColorBrewer")
library("gridExtra")
library("igraph")
library("StAMPP")
library("lattice")
library("reshape2")
library("ggrepel") # for text labels
library("pcadapt")
library("data.table")
library("ggthemes")
library("ggpubr")
library("boot")
library("plotly")
library("forcats")
library("plotly")
Here, we have 43 samples of Varroa mites for which we have sequenced the whole genome succesfully. We use the leaflet package to create an interactive map showing the location and from which host, each Varroa mite was collected.
ID = Name code of the samples obtained from the CSIRO varroa collection
Species = Identified species according to mtDNA of the mtDNA COX1 partial gene
Host collected = Varroa mites were collected from within honey bee colonies identified by the collector as either Western honey bee Apis mellifera or Eastern honey bee Apis cerana.
Location/Country/Date = Exact informations on sampling
Site = Code name of each sampling site corresponding on the map Approximate.X/Y = GPS coordinates from samples obtained before 2008 are infered from approximate location given. All samples geo-coordinates after 2008 were obtained directly from field collection.
COI_haplogroup = mtDNA haplogroup identified from mtDNA COX1 partial region (reconstruct from HiSeq4000 reads and confirmed by Sanger sequencing).
Haplotype = mtDNA haplotype name by concatenating four genes (COX1, COX3, ATP6 and ND2)
Host_read = Host DNA identified by number of reads mapped against either A. cerana or A. mellifera
reads_mellifera/cerana = Number of reads mapped against reference genome of A. cerana or A. mellifera (Q > 20)
# Import the metadata for the Varroa samples
metadata <- read.csv("metadata_varroa_asia_15052020.csv", header = TRUE)
#head(metadata)
## Create a scrolling table
kable(cbind(metadata)) %>%
kable_styling(bootstrap_options = "striped", font_size = 9) %>%
scroll_box(width = "100%", height = "400px") #%>%
| ID | SEQ_ID | Species | Host_collected | Location | Country | Date | Site | Approximate.X | Approximate.Y | COI_haplogroup | Haplotype | Host_read | reads_cerana | reads_mellifera | Mean_Read_Depth | Total_reads | Mapped | Mapping_rate | Mapped_in_pair | Platform |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| V212 | VD212 | V. destructor | A. cerana | Gwangju, Jeonlabug Province | South Korea | 01/1996 | KOR02 | 35.159540 | 126.85260 | VD Korea K1 | VD Korea K1-1/K1-2 | cerana | 56 | 0 | 16.6 | 66272622 | 41762260 | 63.0 | 38388417 | HiSeq 4000 |
| V642_1 | VD642_1 | V. destructor | A. cerana | Zhuhai county, Guangdong Province | China | 11/2001 | CHN08 | 22.614010 | 112.63183 | VD China C1 | VD China C1-3 | cerana | 2857 | 7 | 47.1 | 119107269 | 116682672 | 98.0 | 112230978 | HiSeq 4000 |
| V657_1 | VD657_1 | V. destructor | A. cerana | Zhongshan, Guangdong Province | China | 04/2002 | CHN10 | 23.946090 | 114.69726 | VD China C1 | VD China C1-3 | cerana | 46 | 38 | 21.5 | 54511467 | 53363635 | 97.9 | 49849866 | HiSeq 4000 |
| V651_1 | VD651_1 | V. destructor | A. cerana | Dayao county, Yunnan Province | China | 04/2002 | CHN05 | 25.729510 | 101.33661 | VD China C3 | VD China C3-2 | cerana | 167 | 6 | 28.2 | 71261265 | 69804599 | 98.0 | 66308649 | HiSeq 4000 |
| V577_1 | VD577_1 | V. destructor | A. cerana | Yunnan Agricultural University, Kunming, Yunnan Province | China | 05/2002 | CHN06 | 25.125650 | 102.74318 | VD China C4 | VD China C4-2 | cerana | 125 | 4 | 35.0 | 91721740 | 89577965 | 97.7 | 85757545 | HiSeq 4000 |
| V654_1 | VD654_1 | V. destructor | A. cerana | Xishuangbanna, Yunnan Province | China | 04/2002 | CHN07 | 22.008810 | 100.79714 | VD Viet Nam V1 | VD Viet Nam V1-5 | cerana | 626 | 0 | 25.1 | 63453379 | 62161106 | 98.0 | 58776768 | HiSeq 4000 |
| V658_2 | VD658_2 | V. destructor | A. cerana | Xishuangbanna, Yunnan Province | China | 04/2002 | CHN07 | 22.008810 | 100.79714 | VD Viet Nam V1 | VD Viet Nam V1-6 | cerana | 2693 | 8 | 15.0 | 44339428 | 38654491 | 87.2 | 35187264 | HiSeq 4000 |
| V676_2 | VD676_2 | V. destructor | A. cerana | Chiang Mai University | Thailand | 01/2003 | THA01 | 18.805960 | 98.95336 | VD Viet Nam V1 | VD Viet Nam V1-7 | cerana | 366 | 10 | 5.9 | 29358273 | 19254785 | 65.6 | 17719549 | HiSeq 4000 |
| V149 | VD149 | V. destructor | A. cerana | Phjong, Ninh Binh Province | Viet Nam | 10/1996 | VNM01 | 21.181970 | 106.00161 | VD Viet Nam V1 | VD Viet Nam V1-8 | cerana | 1773 | 76 | 10.7 | 46623619 | 29994005 | 64.3 | 26399304 | HiSeq 4000 |
| V475_1 | VD475_1 | V. destructor | A. cerana | Hoc Chau | Viet Nam | 11/1998 | VNM02 | 20.922080 | 104.75209 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 6 | 555 | 19.3 | 50827244 | 48381793 | 95.2 | 44842843 | HiSeq 4000 |
| V153_2 | VD153_2 | V. destructor | A. mellifera | Jungup | South Korea | 10/1996 | KOR01 | 35.569880 | 126.85589 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 0 | 5824 | 22.2 | 56113114 | 55021332 | 98.1 | 51939126 | HiSeq 4000 |
| V159_1 | VD159_1 | V. destructor | A. mellifera | Jungup | South Korea | 10/1996 | KOR01 | 35.569880 | 126.85589 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 14 | 5707 | 24.2 | 64315412 | 60528716 | 94.1 | 55776173 | HiSeq 4000 |
| V639_1 | VD639_1 | V. destructor | A. mellifera | Heilongjiang Province | China | 09/2001 | CHN01 | 46.617160 | 131.15727 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 20 | 1096 | 17.2 | 45259621 | 43071591 | 95.2 | 39536470 | HiSeq 4000 |
| V625_2 | VD625_2 | V. destructor | A. mellifera | Fengman, Lishu Co., Changchun, Jilin Province | China | 05/2001 | CHN02 | 43.821600 | 126.56227 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 2 | 217 | 16.6 | 46030858 | 41845105 | 90.9 | 37822659 | HiSeq 4000 |
| V622_1 | VD622_1 | V. destructor | A. mellifera | Apiary of Jinzhong, Yuci, Shanxi Province | China | 06/2001 | CHN03 | 37.697790 | 112.70822 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 0 | 1243 | 13.8 | 39024603 | 35070629 | 89.9 | 31675460 | HiSeq 4000 |
| V641_1 | VD641_1 | V. destructor | A. mellifera | Mianzhu County, Sichuan Province | China | 09/2001 | CHN04 | 31.338070 | 104.22074 | VD Korea K1 | VD Korea K1-7 | mellifera | 6 | 871 | 21.0 | 54493356 | 52402973 | 96.2 | 48743391 | HiSeq 4000 |
| V646_1 | VD646_1 | V. destructor | A. mellifera | Zhongshan, Guangdong Province | China | 11/2001 | CHN09 | 23.927990 | 113.15883 | VD Korea K1 | VD Korea K1-7 | mellifera | 4 | 322 | 25.8 | 66768826 | 64203224 | 96.2 | 59914828 | HiSeq 4000 |
| V150_2 | VD150_2 | V. destructor | A. mellifera | Nho Quan | Viet Nam | 10/1996 | VNM03 | 20.322630 | 105.74960 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 26 | 218 | 33.2 | 86978010 | 84333771 | 97.0 | 80000889 | HiSeq 4000 |
| V474_1 | VD474_1 | V. destructor | A. mellifera | Hoc Chau | Viet Nam | 01/1998 | VNM02 | 20.922080 | 104.75209 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 52 | 648 | 24.8 | 63242091 | 61585346 | 97.4 | 58056817 | HiSeq 4000 |
| V028 | VJ028 | V. jacobsoni | A. cerana | UPM, Serdang, Kuala Lumpur | Malaysia | 03/1995 | MYS01 | 2.987890 | 101.73517 | VJ Malaysia 2 | VJ Malaysia 2-1 | cerana | 168 | 66 | 50.4 | 131769709 | 125151456 | 95.0 | 120488777 | HiSeq 4000 |
| V780 | VJ780 | V. jacobsoni | A. cerana | Chauthanh Village, Bentra District | Viet Nam | 03/2003 | VNM06 | 10.304830 | 106.35520 | VJ Laos L1 | VJ Laos L1-3 | cerana | 40372 | 4 | 15.6 | 47957642 | 38953995 | 81.2 | 36880992 | HiSeq 4000 |
| V787_2 | VJ787_2 | V. jacobsoni | A. cerana | Rach Gia Town, Kien Giang | Viet Nam | 08/2004 | VNM04 | 10.021500 | 105.09109 | VJ Laos L1 | VJ Laos L1-4 | cerana | 2325 | 28 | 29.6 | 78330554 | 75988942 | 97.0 | 72622136 | HiSeq 4000 |
| V788_1 (replicate discarded) | VJ788_1 | V. jacobsoni | A. cerana | My tho City, Kien Giang | Viet Nam | 08/2004 | VNM05 | 10.376520 | 106.34388 | VJ Laos L1 | VJ Laos L1-4 | cerana | 1291 | 2 | 7.0 | 43611088 | 19916649 | 45.7 | 18505755 | HiSeq 4000 |
| V788_2 | VJ788_2 | V. jacobsoni | A. cerana | My tho City, Kien Giang | Viet Nam | 08/2004 | VNM05 | 10.376520 | 106.34388 | VJ Laos L1 | VJ Laos L1-4 | cerana | 72 | 2 | 10.8 | 30299193 | 27206968 | 89.8 | 25941038 | HiSeq 4000 |
| V789_1 | VJ789_1 | V. jacobsoni | A. cerana | Giong Trom Town, Ben Tre | Viet Nam | 08/2004 | VNM07 | 10.157180 | 106.50445 | VJ Laos L1 | VJ Laos L1-4 | cerana | 859 | 4 | 23.6 | 71558906 | 58473472 | 81.7 | 56649683 | HiSeq 4000 |
| V333_1 | VJ333_1 | V. jacobsoni | A. cerana | Parung Panjang, Java | Indonesia | 06/1998 | IDN01 | -6.348140 | 106.56935 | VJ Java 1 | VJ Java 1-1 | cerana | 21 | 10 | 17.9 | 46019459 | 44634047 | 97.0 | 42592178 | HiSeq 4000 |
| V341_1 | VJ341_1 | V. jacobsoni | A. cerana | Sukabumi, West Java | Indonesia | 06/1998 | IDN02 | -6.923870 | 106.93156 | VJ Java 1 | VJ Java 1-2 | cerana | 119 | 4 | 24.8 | 63057170 | 61695629 | 97.8 | 59251392 | HiSeq 4000 |
| V347_1 | VJ347_1 | V. jacobsoni | A. cerana | Cililin | Indonesia | 06/1998 | IDN04 | -6.983450 | 107.43608 | VJ Java 1 | VJ Java 1-3 | cerana | 61 | 4 | 18.2 | 48579482 | 45809925 | 94.3 | 42898293 | HiSeq 4000 |
| V363_1 | VJ363_1 | V. jacobsoni | A. cerana | Johor, East Java | Indonesia | 06/1998 | IDN05 | -7.223850 | 112.73320 | VJ Java 1 | VJ Java 1-4 | cerana | 5115 | 64 | 21.3 | 54026529 | 52825280 | 97.8 | 50723792 | HiSeq 4000 |
| V565_2 | VJ565_2 | V. jacobsoni | A. cerana | Gunung Arca, Sukabumi | Indonesia | 04/2002 | IDN03 | -6.926560 | 106.95533 | VJ Java 1 | VJ Java 1-5 | cerana | 714 | 28 | 42.2 | 107959838 | 106090376 | 98.3 | 102687165 | HiSeq 4000 |
| V325_1 | VJ325_1 | V. jacobsoni | A. cerana | Kuta, Lombok | Indonesia | 04/1998 | IDN06 | -8.880270 | 116.28361 | VJ Lombok 2 | VJ Lombok 2-1 | cerana | 152 | 0 | 24.7 | 62599084 | 61346630 | 98.0 | 58922351 | HiSeq 4000 |
| V950_1 | VJ950_1 | V. jacobsoni | A. cerana | Aiyura | Papua New Guinea | 05/1997 | PNG01 | -6.343340 | 145.90809 | VJ PNG 1 | VJ PNG 1-1 | cerana | 475 | 68 | 25.8 | 66335209 | 64182154 | 96.8 | 61215275 | HiSeq 4000 |
| V854_1 | VJ854_1 | V. jacobsoni | A. cerana | Goroka | Papua New Guinea | 05/2008 | PNG04 | -6.083470 | 145.38626 | VJ PNG 1 | VJ PNG 1-1 | cerana | 5230 | 12 | 16.7 | 44613684 | 42039605 | 94.2 | 39498649 | HiSeq 4000 |
| V853_4 | VJ853_4 | V. jacobsoni | A. cerana | Daru | Papua New Guinea | 06/2008 | PNG11 | -9.078190 | 143.20997 | VJ PNG 1 | VJ PNG 1-2 | cerana | 200 | 60 | 10.4 | 35177854 | 32209884 | 91.6 | 29491225 | HiSeq 4000 |
| V983_1 | VJ983_1 | V. jacobsoni | A. cerana | Asaro | Papua New Guinea | 11/2015 | PNG05 | -6.012000 | 145.32500 | VJ PNG 1 | VJ PNG 1-1 | cerana | 284 | 102 | 19.1 | 50283170 | 47881277 | 95.2 | 44989750 | HiSeq 4000 |
| V847 | VJ847 | V. jacobsoni | A. mellifera | Wabag, Enga Province | Papua New Guinea | 06/2008 | PNG10 | -5.482670 | 143.66373 | VJ PNG 1 | VJ PNG 1-1 | mellifera | 23 | 115346 | 22.5 | 64504980 | 56218551 | 87.2 | 52755225 | HiSeq 4000 |
| V852_3 | VJ852_3 | V. jacobsoni | A. mellifera | Goroka, Eastern Highlands Province | Papua New Guinea | 05/2008 | PNG04 | -6.083470 | 145.38626 | VJ PNG 1 | VJ PNG 1-1 | mellifera | 4 | 2288 | 20.9 | 55816067 | 54523224 | 97.7 | 51390551 | HiSeq 4000 |
| V857_2 | VJ857_2 | V. jacobsoni | A. mellifera | Benabena, Eastern Highlands Province | Papua New Guinea | 05/2008 | PNG03 | -6.128125 | 145.42904 | VJ PNG 1 | VJ PNG 1-1 | mellifera | 20 | 429 | 18.8 | 57477057 | 46759609 | 81.4 | 44551404 | HiSeq 4000 |
| V856_1 | VJ856_1 | V. jacobsoni | A. mellifera | Henganofi Border, Eastern Highlands Province | Papua New Guinea | 05/2008 | PNG02 | -6.260350 | 145.59681 | VJ PNG 1 | VJ PNG 1-1 | mellifera | 15 | 9803 | 13.6 | 58525652 | 37045620 | 63.3 | 31803798 | HiSeq 4000 |
| V952_1 | VJ952_1 | V. jacobsoni | A. mellifera | Amayufa, Eastern Highlands Province | Papua New Guinea | 05/2014 | PNG06 | -5.975581 | 145.27676 | VJ PNG 1 | VJ PNG 1-1 | mellifera | 18 | 36 | 22.3 | 61969612 | 55519034 | 89.6 | 52870625 | HiSeq 4000 |
| V953_2 | VJ953_2 | V. jacobsoni | A. mellifera | Benabena, Eastern Highlands Province | Papua New Guinea | 05/2014 | PNG03 | -6.128125 | 145.42904 | VJ PNG 1 | VJ PNG 1-1 | mellifera | 22 | 26196 | 22.3 | 62054570 | 56103596 | 90.4 | 53032044 | HiSeq 4000 |
| V954_2 | VJ954_2 | V. jacobsoni | A. mellifera | Windy Lodge | Papua New Guinea | 06/2014 | PNG07 | -6.035146 | 144.96032 | VJ PNG 1 | VJ PNG 1-1 | mellifera | 14 | 159 | 22.5 | 58553267 | 56180797 | 96.0 | 52846232 | HiSeq 4000 |
| V955_1 | VJ955_1 | V. jacobsoni | A. mellifera | Donakanage | Papua New Guinea | 06/2014 | PNG08 | -5.882240 | 144.82077 | VJ PNG 1 | VJ PNG 1-1 | mellifera | 15 | 460 | 24.7 | 64865555 | 61577522 | 94.9 | 58553870 | HiSeq 4000 |
| V956_4 | VJ956_4 | V. jacobsoni | A. mellifera | Minj | Papua New Guinea | 06/2014 | PNG09 | -5.826640 | 144.40741 | VJ PNG 1 | VJ PNG 1-2 | cerana | 1444 | 69 | 17.1 | 47824500 | 45966716 | 96.1 | 42702569 | HiSeq 4000 |
| Mom #1 | Mom #1 | V. destructor | A. mellifera | OIST Ecology and Evolution experimental apiary, Okinawa | Japan | 02/2018 | JPN01 | 26.456448 | 127.82386 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 0 | 7448 | 20.2 | 52944228 | 50821794 | 96.0 | 52944228 | HiSeq 4000 |
| Son #1 | Son #1 | V. destructor | A. mellifera | OIST Ecology and Evolution experimental apiary, Okinawa | Japan | 02/2018 | JPN01 | 26.456448 | 127.82386 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 0 | 341 | 19.8 | 51620028 | 49266370 | 95.4 | 51620028 | HiSeq 4000 |
| Mom #2 | Mom #2 | V. destructor | A. mellifera | OIST Ecology and Evolution experimental apiary, Okinawa | Japan | 02/2018 | JPN01 | 26.456448 | 127.82386 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 0 | 50765 | 24.8 | 65220912 | 62435130 | 95.7 | 65220912 | HiSeq 4000 |
| Son #2 | Son #2 | V. destructor | A. mellifera | OIST Ecology and Evolution experimental apiary, Okinawa | Japan | 02/2018 | JPN01 | 26.456448 | 127.82386 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 0 | 422 | 18.7 | 48499085 | 46511492 | 95.9 | 48499085 | HiSeq 4000 |
| Mom #7 | Mom #7 | V. destructor | A. mellifera | OIST Ecology and Evolution experimental apiary, Okinawa | Japan | 02/2018 | JPN01 | 26.456448 | 127.82386 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 0 | 563 | 16.3 | 43293860 | 41401282 | 95.6 | 43293860 | HiSeq 4000 |
| Son #7 | Son #7 | V. destructor | A. mellifera | OIST Ecology and Evolution experimental apiary, Okinawa | Japan | 02/2018 | JPN01 | 26.456448 | 127.82386 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 0 | 660 | 18.8 | 49134877 | 46899128 | 95.5 | 49134877 | HiSeq 4000 |
| Mom #15 | Mom #15 | V. destructor | A. mellifera | OIST Ecology and Evolution experimental apiary, Okinawa | Japan | 05/2018 | JPN01 | 26.456448 | 127.82386 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 0 | 100 | 22.5 | 59835146 | 56138955 | 93.8 | 59835146 | HiSeq 4000 |
| Son #15 | Son #15 | V. destructor | A. mellifera | OIST Ecology and Evolution experimental apiary, Okinawa | Japan | 05/2018 | JPN01 | 26.456448 | 127.82386 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 1 | 193 | 21.4 | 57025270 | 53922121 | 94.6 | 57025270 | HiSeq 4000 |
| Mom #17 | Mom #17 | V. destructor | A. mellifera | OIST Ecology and Evolution experimental apiary, Okinawa | Japan | 05/2018 | JPN01 | 26.456448 | 127.82386 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 2 | 40283 | 22.1 | 58014781 | 55352009 | 95.4 | 58014781 | HiSeq 4000 |
| Son #17 | Son #17 | V. destructor | A. mellifera | OIST Ecology and Evolution experimental apiary, Okinawa | Japan | 05/2018 | JPN01 | 26.456448 | 127.82386 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 2 | 949 | 33.0 | 87139528 | 82953771 | 95.2 | 87139528 | HiSeq 4000 |
| Mom #19 | Mom #19 | V. destructor | A. mellifera | OIST Ecology and Evolution experimental apiary, Okinawa | Japan | 05/2018 | JPN01 | 26.456448 | 127.82386 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 0 | 24181 | 18.1 | 47748745 | 45343415 | 95.0 | 47748745 | HiSeq 4000 |
| Son #19 | Son #19 | V. destructor | A. mellifera | OIST Ecology and Evolution experimental apiary, Okinawa | Japan | 05/2018 | JPN01 | 26.456448 | 127.82386 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 5 | 498 | 24.5 | 65304925 | 61543197 | 94.2 | 65304925 | HiSeq 4000 |
#row_spec(1:23, bold = T, color = "white", background = "#E60000") %>%
#row_spec(24:39, bold = T, color = "white", background = "#FF9999") %>%
#row_spec(40:74, bold = T, color = "white", background = "#0000FF") %>%
#row_spec(75:96, bold = T, color = "white", background = "#3399FF")
# Download the country borders layer
#download.file("http://thematicmapping.org/downloads/TM_WORLD_BORDERS_SIMPL-0.3.zip" , destfile="world_shape_file.zip")
#system("unzip world_shape_file.zip")
world_spdf=readOGR(dsn= getwd() , layer="TM_WORLD_BORDERS_SIMPL-0.3")
## OGR data source with driver: ESRI Shapefile
## Source: "/Users/alphamanae/Documents/GitHub/varroa-host-jump/R_data", layer: "TM_WORLD_BORDERS_SIMPL-0.3"
## with 246 features
## It has 11 fields
## Integer64 fields read as strings: POP2005
data.vdac <- metadata %>% filter(Species == "V. destructor") %>% filter(Host_collected == "A. cerana")
data.vdam <- metadata %>% filter(Species == "V. destructor") %>% filter(Host_collected == "A. mellifera")
data.vjac <- metadata %>% filter(Species == "V. jacobsoni") %>% filter(Host_collected == "A. cerana")
data.vjam <- metadata %>% filter(Species == "V. jacobsoni") %>% filter(Host_collected == "A. mellifera")
The following interactive map showed the distribution of both V. destructor and V. jacobsoni samples on either their original host the eastern honey bee A. cerana or on the new host the western honey bee A. mellifera.
If you pass your mouse over a particular point, a pop-up will appear presenting essential information about the particular Varroa mite sample. On the right top corner, a layer control button allow to subselect species/host couple.
V. destructor on A. cerana
V. destructor on A. mellifera
V. jacobsoni on A. cerana
V. jacobsoni on A. mellifera
## Create the interactive map with leaflet
leaflet(metadata) %>%
addTiles(group = "OSM (default)") %>%
## add the layer other than default we would like to use for background
addProviderTiles(providers$CartoDB.PositronNoLabels, group = "Positron NoLabels") %>%
## adding the layers with V. destructor
addCircleMarkers(data = data.vdac, data.vdac$Approximate.Y, data.vdac$Approximate.X,
weight = 0.5,
col = "#FB0000",
radius = 4,
fillOpacity = 0.9,
stroke = T,
popup = ~paste(sep = "<br/>",
paste ("Name: ", as.character(ID)),
paste ("Host: ", as.character(Host_collected)),
paste ("Date: ", as.character(Date)),
paste ("mtDNA: ", as.character(Haplotype))),
group = 'V. destructor on A. cerana') %>%
addCircleMarkers(data = data.vdam, data.vdam$Approximate.Y, data.vdam$Approximate.X,
weight = 0.5,
col = "#FF9999",
radius = 4,
fillOpacity = 0.9,
stroke = T,
popup = ~paste(sep = "<br/>",
paste ("Name: ", as.character(ID)),
paste ("Host: ", as.character(Host_collected)),
paste ("Date: ", as.character(Date)),
paste ("mtDNA: ", as.character(Haplotype))),
group = 'V. destructor on A. mellifera') %>%
## adding the layers with V. jacobsoni
addCircleMarkers(data = data.vjac, data.vjac$Approximate.Y, data.vjac$Approximate.X,
weight = 0.5,
col = "#0000FF",
radius = 4,
fillOpacity = 0.9,
stroke = T,
popup = ~paste(sep = "<br/>",
paste ("Name: ", as.character(ID)),
paste ("Host: ", as.character(Host_collected)),
paste ("Date: ", as.character(Date)),
paste ("mtDNA: ", as.character(Haplotype))),
group = 'V. jacobsoni on A. cerana') %>%
addCircleMarkers(data = data.vjam, data.vjam$Approximate.Y, data.vjam$Approximate.X,
weight = 0.5,
col = "#00CCFF",
radius = 4,
fillOpacity = 0.9,
stroke = T,
popup = ~paste(sep = "<br/>",
paste ("Name: ", as.character(ID)),
paste ("Host: ", as.character(Host_collected)),
paste ("Date: ", as.character(Date)),
paste ("mtDNA: ", as.character(Haplotype))),
group = 'V. jacobsoni on A. mellifera') %>%
## adding the control button to remove or add layers of points
addLayersControl(position = "topright",
baseGroups = c("OSM (default)", "Positron NoLabels"),
overlayGroups = c("V. destructor on A. cerana",
"V. destructor on A. mellifera",
"V. jacobsoni on A. cerana",
"V. jacobsoni on A. mellifera"),
options = layersControlOptions(collapsed = TRUE)) %>%
## show the positron background prerably to the OSM layer
showGroup("Positron NoLabels")
Here we plot the mean read depth computed from each .bam file mapped to the reference vdes3.0.fasta using samtools depth -a FILE.bam | awk '{c++;s+=$3}END{print s/c}' in relation to the mapping percentage.
As you can see by moving your pointer on the graph, the individual V788_1 has been removed from the analysis and replaced by the replicate V788_2 with much better reads statistics.
raichu <- ggplot(metadata, aes(x=Mapping_rate, y = Mean_Read_Depth, col=Species, shape=Host_collected, label = ID))
mycol = c("red2", "blue2")
raichu <- raichu + geom_point(size = 3)
raichu <- raichu + scale_color_manual(values = mycol)
raichu <- raichu + scale_size_manual(values = mycol)
raichu <- raichu + theme_calc()
raichu <- raichu + xlab("Mapping rate to Vdes_3.0 reference genome(%)")
raichu <- raichu + ylab("Mean average read depth")
raichu <- raichu + ylim(0,20)
#raichu
## make an interactive version of the scatter plot
intraichu <- ggplotly(raichu)
intraichu
We highlight the mean read depth regarding the general distribution, where is each group of individuals “species per host” located in it.
evoli <- ggplot(metadata, aes(Mean_Read_Depth, fill = Species))
evoli <- evoli + geom_histogram(bins = 44)
evoli <- evoli + gghighlight()
evoli <- evoli + facet_wrap(~Species~Host_collected)
evoli <- evoli + scale_fill_manual(values = c("red", "blue"))
evoli <- evoli + theme_classic()
evoli <- evoli + labs(x = "Mean Read Depth using NextGenMap", y = "Number of samples")
evoli <- evoli + ggtitle("Read depth does not appear as unbalanced accross species and host")
evoli
Following is a plot of the total number of reads generated per individual and the proportion of reads mapped onto the reference.
togepi <- ggplot(metadata, aes(x=ID, label = Mapping_rate))
togepi <- togepi + geom_bar(aes(weight=Total_reads), fill="black", position="dodge")
togepi <- togepi + geom_bar(aes(weight=Mapped), fill="#1d96f2", position="dodge")
togepi <- togepi + theme(axis.text.x = element_text(angle = 90, hjust = 1))
togepi <- togepi + labs(x = "Individuals", y = "Number of reads")
togepi <- togepi + ggtitle("Proportion of reads mapped against reference using NextGenMap")
togepi
Before proceeding to DNA extraction with a destructive method, we measured the body size of each Varroa mite whenever possible on both dorsal and ventral view. We measured seven morphological characters but retained only the body width and length for comparison here between Varroa species and their associated honey bee hosts.
The scatterplot also shows the average values obtained by Anderson and Trueman (2000) Exp. App. Acar. for V. destructor, V. jacobsoni and undetermined mites from Luzon and Mindanao.
# Table with 512 morpho data for the moment
morpho <- read.csv("Varroa-morpho-subset.csv", header = TRUE)
#kable(cbind(morpho)) %>%
# kable_styling(bootstrap_options = "striped", font_size = 10) %>%
# scroll_box(width = "100%", height = "400px")
morpho %>%
plot_ly(x = ~BW_V, y = ~BL_V,
color = ~Species,
symbol = ~Host,
hoverinfo = "text",
text = ~paste("Name: ", Sample.ID, "<br>",
"Species: ", Species, "<br>",
"Host: ", Host, "<br>",
"Country: ", Country.Island)) %>%
add_markers(colors = c("indianred2", "red3", "skyblue2", "blue2", "black"),
symbols = c( "triangle-up", "circle", "x"),
size = 3) %>%
layout(xaxis = list(title = "Body width (mm)"),
yaxis = list(title = "Body length (mm)"),
title = "Variation in body size measured from ventral view")
morpho %>%
plot_ly(x = ~BW_D, y = ~BL_D,
color = ~Species,
symbol = ~Host,
hoverinfo = "text",
text = ~paste("Name: ", Sample.ID, "<br>",
"Species: ", Species, "<br>",
"Host: ", Host, "<br>",
"Country: ", Country.Island)) %>%
add_markers(colors = c("indianred2", "red3", "skyblue2", "blue2", "black"),
symbols = c( "triangle-up", "circle", "x"),
size = 3) %>%
layout(xaxis = list(title = "Body width (mm)"),
yaxis = list(title = "Body length (mm)"),
title = "Variation in body size measured from dorsal view")
STATS TEST?
After formatting the variant call file .vcf into a plink format by renaming each chromosome (e.g., NW_019211454.1 by 1, NW_019211455.1 by 2, …), and applying the option pca, we obtained the two files .eigenval and .eigenvec.
The following 2D and 3D PCAs are based on the ~1.4 million biallelic SNPs with a minor allelic frequency higher than 5%.
pca.all<- read_table2("43indplink.eigenvec", col_names = FALSE)
eigenval.all <- scan("43indplink.eigenval")
## Here the first two columns contains the individuals names as we did not specify it in plink earlier
# We remove the first column and assign the new name of column 1
pca.all <- pca.all[,-1]
names(pca.all)[1] <- "ID_NAMES"
# Rename the columns for PCA axis
names(pca.all)[2:ncol(pca.all)] <- paste0("PC", 1:(ncol(pca.all)-1))
## Joint both table by using the FILE_ID as common parameters
metapca.all <- left_join(pca.all, metadata, by = c("ID_NAMES" = "SEQ_ID"))
### PLOTS
# first convert to percentage variance explained
pve <- data.frame(PC = 1:20, pve = eigenval.all/sum(eigenval.all)*100)
# make Eigenplot
a <- ggplot(pve, aes(PC, pve)) + geom_bar(stat = "identity")
a + ylab("Percentage variance explained") + theme_light()
metapca.all %>%
plot_ly(x = ~PC1, y = ~PC2,
color = ~Species,
symbol = ~Host_collected,
hoverinfo = "text",
text = ~paste("Name: ", ID, "<br>",
"Species: ", Species, "<br>",
"Host: ", Host_collected, "<br>",
"Country: ", Country)) %>%
add_markers(colors = c("red3", "blue2"),
symbols = c( "triangle-up", "circle"),
size = 4)%>%
layout(title = "Varroa mites genetically differentiate between hosts in both species")
metapca.all %>%
plot_ly(x = ~PC1, y = ~PC2, z = ~PC3,
color = ~Country,
hoverinfo = "text",
text = ~paste("Name: ", ID, "<br>",
"Species: ", Species, "<br>",
"Host: ", Host_collected, "<br>",
"Country: ", Country,
"Haplogroup: ", COI_haplogroup)) %>%
add_markers(colors = c("red3", "purple", "blue2", "green", "orange", "brown", "pink"),
size = 10) %>%
layout(title = "Mites on A. cerana genetically differs within the native range")
metapca.all %>%
plot_ly(x = ~PC1, y = ~PC2, z = ~PC3,
color = ~COI_haplogroup,
hoverinfo = "text",
text = ~paste("Name: ", ID, "<br>",
"Species: ", Species, "<br>",
"Host: ", Host_collected, "<br>",
"Country: ", Country,
"Haplogroup: ", COI_haplogroup)) %>%
add_markers(colors = c("red3", "purple", "blue2", "green", "orange", "brown", "pink"),
size = 10) %>%
layout(title = "Mitochondrial lineages can be a pre-indicator of whole genome differentiation")
Same idea as the previous PCA, except that only 236,713 SNPs were kept to reduce the effect of linkage disequilibrium.
We can see that even with a smaller reduce number of SNPs, Varroa mites species are genetically distinct and
pca.allLD<- read_table2("43ind-ldpruned.eigenvec", col_names = FALSE)
eigenval.allLD <- scan("43ind-ldpruned.eigenval")
## Here the first two columns contains the individuals names as we did not specify it in plink earlier
# We remove the first column and assign the new name of column 1
pca.allLD <- pca.allLD[,-1]
names(pca.allLD)[1] <- "ID_NAMES"
# Rename the columns for PCA axis
names(pca.allLD)[2:ncol(pca.allLD)] <- paste0("PC", 1:(ncol(pca.allLD)-1))
## Joint both table by using the FILE_ID as common parameters
metapca.allLD <- left_join(pca.allLD, metadata, by = c("ID_NAMES" = "SEQ_ID"))
### PLOTS
# first convert to percentage variance explained
pve <- data.frame(PC = 1:20, pve = eigenval.allLD/sum(eigenval.allLD)*100)
# make Eigenplot
a <- ggplot(pve, aes(PC, pve)) + geom_bar(stat = "identity")
a + ylab("Percentage variance explained") + theme_light()
metapca.allLD %>%
plot_ly(x = ~PC1, y = ~PC2,
color = ~Species,
symbol = ~Host_collected,
hoverinfo = "text",
text = ~paste("Name: ", ID, "<br>",
"Species: ", Species, "<br>",
"Host: ", Host_collected, "<br>",
"Country: ", Country)) %>%
add_markers(colors = c("red3", "blue2"),
symbols = c( "triangle-up", "circle"),
size = 4) %>%
layout(title = "~200,000 LD pruned SNPs still allowed to see host genetic divergence")
metapca.allLD %>%
plot_ly(x = ~PC1, y = ~PC2, z = ~PC3,
color = ~Country,
hoverinfo = "text",
text = ~paste("Name: ", ID, "<br>",
"Species: ", Species, "<br>",
"Host: ", Host_collected, "<br>",
"Country: ", Country,
"Haplogroup: ", COI_haplogroup)) %>%
add_markers(colors = c("red3", "purple", "blue2", "green", "orange", "brown", "pink"),
size = 10) %>%
layout(title = "3D PCA color-coded by country")
metapca.allLD %>%
plot_ly(x = ~PC1, y = ~PC2, z = ~PC3,
color = ~COI_haplogroup,
hoverinfo = "text",
text = ~paste("Name: ", ID, "<br>",
"Species: ", Species, "<br>",
"Host: ", Host_collected, "<br>",
"Country: ", Country,
"Haplogroup: ", COI_haplogroup)) %>%
add_markers(colors = c("red3", "purple", "blue2", "green", "orange", "brown", "pink"),
size = 10) %>%
layout(title = "3D PCA color-coded by mtDNA COX1 haplogroup")
After estimating the demographic history for each Varroa species on their novel and original hosts using fastsimcoal2, we plotted here the distribution of the log likelihood from the 100 run replicates for each scenario and each SFS subset size.
This method for model selection comes on top of the AIC calculation and is followed the same methods as described in Meier et al. 2016. Mol. Ecol.. For the AIC values, see Table S7.
vd.lhoods1050 <- read.csv("fsc_run/lhoods_1050SNPS_vdno475.csv", header = TRUE)
vd.lhoods10500 <- read.csv("fsc_run/lhoods_10500SNPS_vdno475.csv", header = TRUE)
vd.lhoods35000 <- read.csv("fsc_run/lhoods_35000SNPS_vdno475.csv", header = TRUE)
vd.lhoods105000 <- read.csv("fsc_run/lhoods_105000SNPS_vdno475.csv", header = TRUE)
par(mfrow=c(1,4))
#Plot 1 for likelihoods with smaller subset
boxplot(range = 0,vd.lhoods1050$scenario1_div [3:102], vd.lhoods1050$scenario2_divGwt2n[3:102], vd.lhoods1050$scenario3_divGwt[3:102], vd.lhoods1050$scenario4_IM[3:102], vd.lhoods1050$scenario5_IMGwt2n[3:102], vd.lhoods1050$scenario6_IMGwt[3:102], ylab="log 10 Likelihood",xaxt="n")
points(1,vd.lhoods1050$scenario1_div [2],col = "red", pch = 20)
points(2,vd.lhoods1050$scenario2_divGwt2n [2],col = "red", pch = 20)
points(3,vd.lhoods1050$scenario3_divGwt [2],col = "red", pch = 20)
points(4,vd.lhoods1050$scenario4_IM [2],col = "red", pch = 20)
points(5,vd.lhoods1050$scenario5_IMGwt2n [2],col = "red", pch = 20)
points(6,vd.lhoods1050$scenario6_IMGwt [2],col = "red", pch = 20)
axis(side=1,at=1:6, labels=c("DIV1","DIV2","DIV3","IM4","IM5","IM6"))
title("V. destructor 377 SNPS")
#Plot 2 for likelihoods with second smaller subset
boxplot(range = 0,vd.lhoods10500$scenario1_div [3:102], vd.lhoods10500$scenario2_divGwt2n[3:102], vd.lhoods10500$scenario3_divGwt[3:102], vd.lhoods10500$scenario4_IM[3:102], vd.lhoods10500$scenario5_IMGwt2n[3:102], vd.lhoods10500$scenario6_IMGwt[3:102], ylab="log 10 Likelihood",xaxt="n")
points(1,vd.lhoods10500$scenario1_div [2],col = "red", pch = 20)
points(2,vd.lhoods10500$scenario2_divGwt2n [2],col = "red", pch = 20)
points(3,vd.lhoods10500$scenario3_divGwt [2],col = "red", pch = 20)
points(4,vd.lhoods10500$scenario4_IM [2],col = "red", pch = 20)
points(5,vd.lhoods10500$scenario5_IMGwt2n [2],col = "red", pch = 20)
points(6,vd.lhoods10500$scenario6_IMGwt [2],col = "red", pch = 20)
axis(side=1,at=1:6, labels=c("DIV1","DIV2","DIV3","IM4","IM5","IM6"))
title("V. destructor 3,760 SNPS")
#Plot 3 for likelihoods with third subset
boxplot(range = 0,vd.lhoods35000$scenario1_div [3:102], vd.lhoods35000$scenario2_divGwt2n[3:102], vd.lhoods35000$scenario3_divGwt[3:102], vd.lhoods35000$scenario4_IM[3:102], vd.lhoods35000$scenario5_IMGwt2n[3:102], vd.lhoods35000$scenario6_IMGwt[3:102], ylab="log 10 Likelihood",xaxt="n")
points(1,vd.lhoods35000$scenario1_div [2],col = "red", pch = 20)
points(2,vd.lhoods35000$scenario2_divGwt2n [2],col = "red", pch = 20)
points(3,vd.lhoods35000$scenario3_divGwt [2],col = "red", pch = 20)
points(4,vd.lhoods35000$scenario4_IM [2],col = "red", pch = 20)
points(5,vd.lhoods35000$scenario5_IMGwt2n [2],col = "red", pch = 20)
points(6,vd.lhoods35000$scenario6_IMGwt [2],col = "red", pch = 20)
axis(side=1,at=1:6, labels=c("DIV1","DIV2","DIV3","IM4","IM5","IM6"))
title("V. destructor 12,665 SNPS")
#Plot 4 for likelihoods with larger subset
boxplot(range = 0,vd.lhoods105000$scenario1_div [3:102], vd.lhoods105000$scenario2_divGwt2n[3:102], vd.lhoods105000$scenario3_divGwt[3:102], vd.lhoods105000$scenario4_IM[3:102], vd.lhoods105000$scenario5_IMGwt2n[3:102], vd.lhoods105000$scenario6_IMGwt[3:102], ylab="log 10 Likelihood",xaxt="n")
points(1,vd.lhoods105000$scenario1_div [2],col = "red", pch = 20)
points(2,vd.lhoods105000$scenario2_divGwt2n [2],col = "red", pch = 20)
points(3,vd.lhoods105000$scenario3_divGwt [2],col = "red", pch = 20)
points(4,vd.lhoods105000$scenario4_IM [2],col = "red", pch = 20)
points(5,vd.lhoods105000$scenario5_IMGwt2n [2],col = "red", pch = 20)
points(6,vd.lhoods105000$scenario6_IMGwt [2],col = "red", pch = 20)
axis(side=1,at=1:6, labels=c("DIV1","DIV2","DIV3","IM4","IM5","IM6"))
title("V. destructor 38,108 SNPS")
After summarizing the 100 best replicates from the 100 simulated SFS run, we computed the mean estimate value and 95% confidence interval for each demographic parameter using boot package.
my.mean = function(x, indices) {
return( mean( x[indices] ) )
}
# import the file containing the 100 best runs
vd.boot6 <- read.csv("fsc_run/vdno475_105000nopv_BOOT.csv", header = TRUE)
# Effective size for Varroa mites on original host Apis cerana
N_VdAc1.boot = boot(vd.boot6$NVDAC1, my.mean, 10000)
my.mean(vd.boot6$NVDAC1,1:length(vd.boot6$NVDAC1))
## [1] 412.5
boot.ci(N_VdAc1.boot)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = N_VdAc1.boot)
##
## Intervals :
## Level Normal Basic
## 95% (400.6, 424.4 ) (400.7, 424.1 )
##
## Level Percentile BCa
## 95% (400.9, 424.3 ) (400.8, 424.3 )
## Calculations and Intervals on Original Scale
# Effective size for Varroa mites on new host Apis mellifera
N_VdAm0.boot = boot(vd.boot6$NVDAM0, my.mean, 10000)
my.mean(vd.boot6$NVDAM0,1:length(vd.boot6$NVDAM0))
## [1] 568.11
boot.ci(N_VdAm0.boot)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = N_VdAm0.boot)
##
## Intervals :
## Level Normal Basic
## 95% (550.2, 586.2 ) (550.4, 586.6 )
##
## Level Percentile BCa
## 95% (549.6, 585.8 ) (549.5, 585.7 )
## Calculations and Intervals on Original Scale
# Number of generations when new population of Varroa mites was founded
TBOT.boot = boot(vd.boot6$TBOT, my.mean, 10000)
my.mean(vd.boot6$TBOT,1:length(vd.boot6$TBOT))
## [1] 722.98
boot.ci(TBOT.boot)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = TBOT.boot)
##
## Intervals :
## Level Normal Basic
## 95% (690.0, 755.8 ) (689.9, 755.7 )
##
## Level Percentile BCa
## 95% (690.3, 756.1 ) (690.4, 756.1 )
## Calculations and Intervals on Original Scale
# Growth rate of Varroa mites on new host Apis mellifera
GAM.boot = boot(vd.boot6$GAM, my.mean, 10000)
my.mean(vd.boot6$GAM,1:length(vd.boot6$GAM))
## [1] -0.0034098
boot.ci(GAM.boot)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = GAM.boot)
##
## Intervals :
## Level Normal Basic
## 95% (-0.0036, -0.0032 ) (-0.0036, -0.0032 )
##
## Level Percentile BCa
## 95% (-0.0036, -0.0032 ) (-0.0036, -0.0032 )
## Calculations and Intervals on Original Scale
# Effective size for Varroa mites able to jump on A. mellifera to create new population
JUMP.boot = boot(vd.boot6$JUMP, my.mean, 10000)
my.mean(vd.boot6$JUMP,1:length(vd.boot6$JUMP))
## [1] 165.04
boot.ci(JUMP.boot)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = JUMP.boot)
##
## Intervals :
## Level Normal Basic
## 95% (160.3, 169.7 ) (160.3, 169.6 )
##
## Level Percentile BCa
## 95% (160.5, 169.8 ) (160.5, 169.8 )
## Calculations and Intervals on Original Scale
N0M1.boot = boot(vd.boot6$N0M1, my.mean, 10000)
my.mean(vd.boot6$N0M1,1:length(vd.boot6$N0M1))
## [1] 0.0005941365
boot.ci(N0M1.boot)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = N0M1.boot)
##
## Intervals :
## Level Normal Basic
## 95% ( 0.0003, 0.0009 ) ( 0.0003, 0.0008 )
##
## Level Percentile BCa
## 95% ( 0.0003, 0.0009 ) ( 0.0004, 0.0011 )
## Calculations and Intervals on Original Scale
N1M0.boot = boot(vd.boot6$N1M0, my.mean, 10000)
my.mean(vd.boot6$N1M0,1:length(vd.boot6$N1M0))
## [1] 0.7604321
boot.ci(N1M0.boot)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = N1M0.boot)
##
## Intervals :
## Level Normal Basic
## 95% ( 0.7586, 0.7623 ) ( 0.7586, 0.7623 )
##
## Level Percentile BCa
## 95% ( 0.7586, 0.7623 ) ( 0.7586, 0.7623 )
## Calculations and Intervals on Original Scale
MACtoAM.boot = boot(vd.boot6$MACtoAM, my.mean, 10000)
my.mean(vd.boot6$MACtoAM,1:length(vd.boot6$MACtoAM))
## [1] 0.001885131
boot.ci(MACtoAM.boot)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = MACtoAM.boot)
##
## Intervals :
## Level Normal Basic
## 95% ( 0.0018, 0.0019 ) ( 0.0018, 0.0019 )
##
## Level Percentile BCa
## 95% ( 0.0018, 0.0019 ) ( 0.0018, 0.0019 )
## Calculations and Intervals on Original Scale
MAMtoAC.boot = boot(vd.boot6$MAMtoAC, my.mean, 10000)
my.mean(vd.boot6$MAMtoAC,1:length(vd.boot6$MAMtoAC))
## [1] 1.044978e-06
boot.ci(MAMtoAC.boot)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = MAMtoAC.boot)
##
## Intervals :
## Level Normal Basic
## 95% ( 0, 0 ) ( 0, 0 )
##
## Level Percentile BCa
## 95% ( 0, 0 ) ( 0, 0 )
## Calculations and Intervals on Original Scale
Here we modified the .par file from the best scenario 6 for V. destructor by adding the mean values computed in the step before. Using the ParFileInterpreter-v6.3.1.r developed by Victor Sousa and available on fastsimcoal2 website, we can now visulized the best demographic scenario.
#...........................................................................................
# (c) Laurent Excoffier and Vitor Sousa June-November 2015
#
# Small R program to draw the evolutonary scenario described by a given par file
# This is mainly for visual checking that the modeled scenario corresponds to
# what was intended
#
#...........................................................................................
#
args=commandArgs(TRUE)
print(args)
## character(0)
#Expects par file name on command line
if(length(args)) {
parFileName=args[1]
} else {
#REPLACE BY THE NAME OF THE PAR FILE YOU WANT TO ANALYSE
parFileName="fsc_run/VDM6_downsample_IMGwt_49_boot.par"
}
separator=" "
migrMatCol="coral"
admixCol="blue"
popFusionColor="black"
popCol="lightgrey"
popBorderCol="black"
inbreedColor="lightblue"
oldSampColor="olivedrab4"
timeCol="tan4"
growthCol="hotpink3"
propLastsegment=0.05
migMatNameProp=0.8
migMatLineLength=0.3
timeProp=0.6
maxRadius=1/40
minRadius=maxRadius/3
arrowLength=0.2
logScaleAxis=""
timeOffset=0.25
migrOffset=0.1
curvedArrowLTY=1
drawLogPopSize=T
plotMigrRates=T
migrRateTextSizeCEX=0.5
#Define plot area size for PDF
pdf.x.size=7
pdf.y.size=10
rescalingFactor=1.0 #Don't touch this!
printPDF=F
# parFile=readLines(con=parFileName)
########################### Reading par file #########################
parFile=scan(parFileName, character(0), sep = "\n", strip.white = TRUE) # separate each line
#pdfFileName=paste(parFileName, ".pdf", sep="")
#--- Clean input file by removing consecutive separators, and keep ---------------
#--- Function to remove separators within a string
removeTrailingSep=function(string, sep) {
temp=strsplit(string, split=sep)
temp2=temp[[1]][nchar(temp[[1]])>0]
cleanStr=temp2[1]
if (length(temp2)>1) {
for (i in 2:length(temp2)) {
cleanStr=paste(cleanStr, temp2[i], sep=sep)
}
}
cleanStr
}
#--- Replace Keep by -9999
replaceKeep=function(string) {
if (grepl("keep", string)) {
return(gsub("keep", "-9999", string))
}
return(string)
}
#Remove both multiple consecutive whitespace and tabs and replace keep keyword
for (i in 1:length(parFile)) {
parFile[i]=removeTrailingSep(parFile[i], sep='\t')
parFile[i]=removeTrailingSep(parFile[i], sep=' ')
parFile[i]=replaceKeep(parFile[i])
}
#-------------------------------------------------------------------------------
#--- Get number of samples on line 2 -----
l.numsamples=parFile[2]
sp.l=unlist(strsplit(l.numsamples, split=' '))
tab.l=unlist(strsplit(l.numsamples, split='\t'))
if (length(sp.l)>=length(tab.l)) {
numSamples=as.numeric(sp.l[1])
separator=" "
} else {
numSamples=as.numeric(tab.l[1])
separator="\t"
}
#--- Reading numbers on separate lines -----
getNumbers=function(start, parFile, numSamples) {
for (i in 1:numSamples) {
curnum=as.numeric(unlist(strsplit(parFile[start+i], split=separator))[1])
if (i==1) {
num=curnum
} else {
num=c(num, curnum)
}
}
num
}
#--- Get population sizes -----------------
start=3
popSizes=getNumbers(start, parFile, numSamples)
#Rescaling pop sizes
popSizes=round(popSizes*rescalingFactor, digits=0)
iniPopSizes=popSizes
#--- Get sample sizes -----------------
readSampleSizesTimesAndInbreedingLevel=function(start, parFile, numSamples) {
for (i in 1:numSamples) {
curLine=unlist(strsplit(parFile[start+i], split=separator))
curSampSize=as.numeric(curLine[1])
curSampTime=0
curInbreeding=0
if (length(curLine)>1) curSampTime=as.numeric(curLine[2])
if (length(curLine)>2) curInbreeding=as.numeric(curLine[3])
if (i==1) {
sampSize=curSampSize
sampTime=curSampTime
inbreeding=curInbreeding
} else {
sampSize=c(sampSize,curSampSize)
sampTime=c(sampTime,curSampTime)
inbreeding=c(inbreeding,curInbreeding)
}
}
list(ss=sampSize, st=sampTime, inb=inbreeding)
}
start=start+numSamples+1
# sampSizes=getNumbers(start, parFile, numSamples)
sampSizesStats=readSampleSizesTimesAndInbreedingLevel(start, parFile, numSamples)
#Rescaling sample times
sampSizesStats$st=round(sampSizesStats$st*rescalingFactor, digits=0)
sampSizes=sampSizesStats$ss
sampTimes=sampSizesStats$st
inbrCoeff=sampSizesStats$inb
#--- Get growth rates -----------------
start=start+numSamples+1
growthRatesInitial=getNumbers(start, parFile, numSamples)
# save this into growthRates which will be used and updated when printing historical events
growthRates=growthRatesInitial
#--- Get number of migration matrices -----------------
start=start+numSamples+1
numMigMat=as.numeric(unlist(strsplit(parFile[start+1], split=separator))[1])
#--- Read migration matrix
readMigMat=function(start, parFile, numSamples) {
for (i in 1:numSamples) {
curmigs=as.numeric(unlist(strsplit(parFile[start+i], split=separator)))
if (i==1) {
migs=curmigs
} else {
migs=rbind(migs, curmigs)
}
}
rownames(migs)=1:numSamples
migs
}
#--- Get migration matrices as a list --------------
start=start+2
migMats=list()
if (numMigMat>0) {
for (i in 1:numMigMat) {
curMigMat=readMigMat(start, parFile, numSamples)
migMats[[i]]=curMigMat
start=start+numSamples+1
}
}
#Rescaling migration rates
if (numMigMat>0) {
for (i in 1:numMigMat) {
migMats[[i]]=migMats[[i]]/rescalingFactor
}
}
#--- Get number of historical events
start=start+1
numHistEvents=as.numeric(unlist(strsplit(parFile[start], split=separator))[1])
###################### HISTORICAL EVENTS HANDLING ##############################
#..... Read Historical Event .......
last.he.time=0
if (numHistEvents>0) {
for (i in 1:numHistEvents) {
start=start+1
#Take care of nomig keyword
nomig=F
if (grepl("nomig", parFile[start])) {
nomig=T
gsub("nomig", "", parFile[start])
}
curHE=as.numeric(unlist(strsplit(parFile[start], split=separator)))
if (nomig) curHE[7]=-1
#Rescaling time of event
curHE[1]=round(curHE[1]*rescalingFactor, digits=0)
if (i==1) {
histEvents=curHE
last.he.time=curHE[1]
} else {
histEvents=rbind(histEvents, curHE)
if (histEvents[i,1] > last.he.time) last.he.time=histEvents[i,1]
}
}
if (numHistEvents>1) rownames(histEvents)=1:numHistEvents
}
names(last.he.time)=""
last.he.time=as.numeric(last.he.time)
yTimeLimit=0
if (last.he.time!=0) yTimeLimit=last.he.time*(1+propLastsegment) #Add propLastsegment to y axis after last event (to draw stuff)
#Reorder events by their times
if (numHistEvents>1) histEvents=histEvents[order(histEvents[,1],decreasing=FALSE),] else {
if (numHistEvents==1) histEvents=matrix(histEvents, nrow=1, byrow = T)
}
endReadParFile=start
############################## PLOTTING THE MODELED SCENARIO ######################################
#--- Graphical functions ...........................................................................
fullHeadArrow=function(x0, y0, x1, y1, length, angle, color="black", weight=1) {
arrows(x0, y0, x1, y1, length, angle, code=2, lty=1, col=color, lwd=weight)
arrows(x0, y0, x1, y1, length, angle*0.80, code=2, lty=1, col=color, lwd=weight)
arrows(x0, y0, x1, y1, length, angle*0.60, code=2, lty=1, col=color, lwd=weight)
arrows(x0, y0, x1, y1, length, angle*0.40, code=2, lty=1, col=color, lwd=weight)
arrows(x0, y0, x1, y1, length, angle*0.20, code=2, lty=1, col=color, lwd=weight)
arrows(x0, y0, x1, y1, length, angle*0.10, code=2, lty=1, col=color, lwd=weight)
}
drawTriangle=function(growth, x, y, size, aspRatio, color) {
if (growth>0) {
x0=x; y0=y
x1=x-size/2; y1=y+size/2*aspRatio
x2=x+size/2; y2=y+size/2*aspRatio
} else {
x0=x-size/2; y0=y
x1=x; y1=y+size/2*aspRatio
x2=x+size/2; y2=y
}
polygon(c(x0, x1, x2, x0), c(y0, y1, y2, y0), col=color)
return(y+size/2*aspRatio)
}
#--- Computing maximum current pop size ............................................................
maxPopSize=popSizes[1]
minPopSize=popSizes[1]
if (length(popSizes)>1){
for (i in 2:length(popSizes)) {
if (popSizes[i]>maxPopSize) maxPopSize=popSizes[i];
if (popSizes[i]<minPopSize) minPopSize=popSizes[i];
}
}
#--- Find min and max pop sizes over the whole population history ..................................
ps=popSizes
isGrowth=FALSE
if (numHistEvents) {
#Need to keep track ofgrowth rates over time
gRates=growthRates
prevTime=0
for (i in 1:numHistEvents) {
he=histEvents[i,]
curTime=he[1]
sink=he[3]+1
resize=he[5]
growth=he[6]
#Begin by resizing pop sizes due to growth since last event
for (j in 1:length(ps)) {
ps[j]=ps[j]*exp(gRates[j]*(curTime-prevTime))
if (ps[j]>maxPopSize) maxPopSize=ps[j]
if (ps[j]<minPopSize) minPopSize=ps[j]
if (ps[j]==0) print(paste("Deme ", j-1, " reaches size zero at time ", curTime, " due to large negative growth (", gRates[j], ")", sep=""))
if (is.infinite(ps[j])) print(paste("Deme ", j-1, " reaches infinite size at time ", curTime, " due to large positive growth (", gRates[j], ")", sep=""))
}
prevTime=curTime
#Handle transformed keep keyword
if(growth!=-9999) {
gRates[sink]=growth
}
#Implement resize
ps[sink]=ps[sink]*resize
if (ps[sink]>maxPopSize) maxPopSize=ps[sink]
if (ps[sink]<minPopSize) minPopSize=ps[sink]
}
}
#-- Function to compute the circle radius for a given pop size ....................................
interpolRadius=function(curSize, minSize, maxSize, minRadius, maxRadius, logScale) {
if(logScale) {
minSize=log10(minSize)
maxSize=log10(maxSize)
curSize=log10(curSize)
}
curRadius=minRadius+(curSize-minSize)*(maxRadius-minRadius)/(maxSize-minSize)
curRadius
}
#============================== ========================
#============================== BEGIN MAIN PLOT ========================
#============================== ========================
library("plotrix")
library("diagram")
## Loading required package: shape
##
## Attaching package: 'diagram'
## The following object is masked from 'package:sp':
##
## coordinates
if (printPDF) {
pdf(pdfFileName, width=pdf.x.size, height=pdf.y.size)
}
par(xpd=F, mar=c(6,4,3,0.5))
maxRadius=maxRadius*(numSamples+2)
minRadius=minRadius*(numSamples+2)
title=parFileName
plot(x=1:numSamples, type="n", xlab="", ylab="time (gen)", xlim=c(-0.5, numSamples+1.5), cex.main=0.8,
ylim=c(0, yTimeLimit), main=title, xaxt = 'n', log=logScaleAxis, cex.axis=0.9, cex.lab=0.9)
axis(side=1, labels=c("Mig Mat",0:(numSamples-1), " \nTimes"), at=0:(numSamples+1), cex.axis=0.8)
mtext("demes", side=1,line=2, cex=0.8)
w <- par("pin")[1]/diff(par("usr")[1:2])
h <- par("pin")[2]/diff(par("usr")[3:4])
aspRatio <- w/h
# aspRatio=yTimeLimit/(numSamples+2)
#--- Plotting population circles according to their current size and sampling time
for (i in 1:numSamples) {
curRadius=interpolRadius(popSizes[i], minPopSize, maxPopSize, minRadius, maxRadius, drawLogPopSize)
# print(curRadius)
if (i==1) topOfCircle=sampTimes[i]+curRadius*aspRatio else topOfCircle=c(topOfCircle, sampTimes[i]+curRadius*aspRatio)
curColor=popCol
if (sampSizes[i]==0) curColor="white"
if (inbrCoeff[i]==0) {
draw.circle(i, sampTimes[i], radius = curRadius, col=curColor, border=popBorderCol)
} else {
floating.pie(i, sampTimes[i], c(inbrCoeff[i], 1-inbrCoeff[i]), radius = curRadius,
col=c(inbreedColor,curColor), startpos=pi, border=popBorderCol)
}
if (sampTimes[i]>0) {
text(i, sampTimes[i]-curRadius*aspRatio, labels=sampTimes[i], cex=timeProp, col=oldSampColor, pos=1)
}
#Draw a vertical arrow in case of pop growth
curGrowthRate=growthRates[i]
#--- Draw gtrowing or shrinking triangle on top of pop circle to reflect growth type
if (curGrowthRate!=0) {
arLength=0.15
topOfCircle[i]=drawTriangle(curGrowthRate, x=i, y=topOfCircle[i], arLength, aspRatio, color=growthCol)
}
}
#--- Handle first migration matrix .......................
curMigMatNum=0
curvature=0.0075*last.he.time
text(0-migrOffset, 0, labels=0, cex=migMatNameProp, col=migrMatCol)
if (numMigMat) {
curMigMat=migMats[1][[1]]
for (sink in 1:numSamples) {
for (sourc in 1:numSamples) {
if (sink!=sourc & curMigMat[sourc, sink]>0) {
differ=sourc-sink
curvedarrow(from=c(sourc, 0), to=c(sink,0), curve=curvature*(abs(differ)*0.55^abs(differ)), arr.adj=1,
arr.pos=0.5, arr.type="triangle", arr.col=migrMatCol, lwd=1, lty=curvedArrowLTY,
lcol=migrMatCol, arr.length=arrowLength)
if (plotMigrRates) {
curNm=round(curMigMat[sourc, sink]*popSizes[sourc], digits=2)
if (differ>0) {
text(sink+abs(differ)/2, aspRatio*0.15*abs(differ), labels=curNm, cex=migrRateTextSizeCEX, col=migrMatCol)
} else {
text(sourc+abs(differ)/2, -aspRatio*0.15*abs(differ), labels=curNm, cex=migrRateTextSizeCEX, col=migrMatCol)
}
}
}
}
}
}
#---Draw all events on the population tree
lastTime=0
activePops=1:numSamples
numActivePops=numSamples
lastSink=-1
if (numHistEvents) {
for (i in 1:numHistEvents) {
#Extract historical event
he=histEvents[i,]
he.time=he[1]
he.source=he[2]+1 #+1 is due to the use of base 0 for deme number in fsc
he.sink=he[3]+1
he.migr=he[4]
he.resize=he[5]
he.growth=he[6]
if(he.growth==-9999) { #handle transformed keep keyword
he.growth=growthRates[he.sink] #Keep current growth rate
}
he.migrMat=he[7]
if(he.migrMat==-9999) { #handle transformed keep keyword
he.migrMat=curMigMatNum
}
#-- Draw event time
if (i%%2) {
slide=timeOffset
}
else {
slide=-timeOffset
}
#Draw all vertical segments ....................
for (p in 1:length(activePops)) {
#Handle population with older sampling times
if (sampTimes[activePops[p]]>0) {
minTime=max(lastTime, topOfCircle[activePops[p]]);
} else {
minTime=lastTime;
}
if (he.time>sampTimes[activePops[p]] & he.time>topOfCircle[activePops[p]]) {
if (activePops[p]==lastSink | i==1) {
segments(activePops[p], topOfCircle[activePops[p]], activePops[p], he.time)
} else {
segments(activePops[p], minTime, activePops[p], he.time)
}
}
}
#Handle growth rate changes since last event ......................
#Update pop sizes according to current growth rates
if (i>1) {
prev.he=histEvents[i-1,]
branchLength=he.time-lastTime
for (p in 1:length(activePops)) {
curPop=activePops[p]
popSizes[curPop]=popSizes[curPop]*exp(growthRates[curPop]*branchLength);
}
}
#Update growth rate
growthRates[he.sink]=he.growth
lastTime=he.time
#Handle resize of sink pop ........................
popSizes[he.sink]=popSizes[he.sink]*he.resize
curRadius=interpolRadius(popSizes[he.sink], minPopSize, maxPopSize, minRadius, maxRadius, drawLogPopSize)
# popRadius[he.sink]=curRadius
topOfCircle[he.sink]=he.time+curRadius*aspRatio
draw.circle(he.sink, he.time, radius = curRadius, col=popCol, border=popBorderCol)
#--- Draw growing or shrinking triangle on top of pop circle to reflect growth type
if (he.growth!=0) {
arLength=0.15
topOfCircle[he.sink]=drawTriangle(he.growth, x=he.sink, y=topOfCircle[he.sink], arLength, aspRatio, color=growthCol)
}
#Handle population fusion .........................
if (he.migr>=1 & he.sink!=he.source) { #This is a population fusion
if (numActivePops==numSamples) removedPops=he.source else removedPops=c(removedPops, he.source)
numActivePops=numActivePops-1
activePops=(1:numSamples)[-removedPops]
#Draw connecting arrows from source to sink
fullHeadArrow(he.source, he.time, he.sink, he.time, length=0.15, angle=20)
#Redraw time with the right color
text(numSamples+1+slide, he.time, labels=he.time, cex=timeProp, col=popFusionColor)
} else {
#--- Handle admixture event ........................
if (he.migr>0 & he.migr<1) {
#Draw connecting arrows from source to sink
segments(he.source, he.time, he.sink, he.time, col=admixCol, lty=2)
if (he.sink>he.source) {
fullHeadArrow(he.sink-0.15, he.time, he.sink, he.time, length=0.15, angle=20, color=admixCol)
} else {
fullHeadArrow(he.sink+0.15, he.time, he.sink, he.time, length=0.15, angle=20, color=admixCol)
}
#Redraw time with the right color
text(numSamples+1+slide, he.time, labels=he.time, cex=timeProp, col=admixCol)
}
else text(numSamples+1+slide, he.time, labels=he.time, cex=timeProp, col=timeCol)
}
#--- Handle migmat change ............................
if (i!=numHistEvents) nextTime=histEvents[i+1,1] else {
nextTime=yTimeLimit
}
time2DrawArrows=(he.time+nextTime)/2
if (he.migrMat!=curMigMatNum) {
if (he.migrMat>-1) {
curMigMat=migMats[he.migrMat+1][[1]]
for (sink in 1:numSamples) {
for (sourc in 1:numSamples) {
if (sink!=sourc & curMigMat[sourc, sink]>0) {
differ=sourc-sink
curvedarrow(from=c(sourc, time2DrawArrows), to=c(sink, time2DrawArrows), curve=curvature*(abs(differ)*0.55^abs(differ)),
arr.adj=1, arr.pos=0.5, arr.type="triangle", arr.col=migrMatCol, lwd=1, lty=curvedArrowLTY,
lcol=migrMatCol, arr.length=arrowLength)
if (plotMigrRates) {
#Write Nm values
curNm=round(curMigMat[sourc, sink]*popSizes[sourc], digits=2)
if (differ>0) {
text(sink+abs(differ)/2, time2DrawArrows+aspRatio*0.15*abs(differ), labels=curNm, cex=migrRateTextSizeCEX, col=migrMatCol)
} else {
text(sourc+abs(differ)/2, time2DrawArrows-aspRatio*0.15*abs(differ), labels=curNm, cex=migrRateTextSizeCEX, col=migrMatCol)
}
}
}
}
# curMigMatNum=he.migrMat
}
}
#--- Draw separation between migration matrices numbers on the left
segments(-migMatLineLength/2, he.time, migMatLineLength/2, he.time, lty=3, col=migrMatCol)
}
curMigMatNum=he.migrMat
#Output current valid migration matrix
if (i%%2) {
slide=migrOffset
}
else {
slide=-migrOffset
}
if (he.migrMat>-1) {
migText=he.migrMat
curCex=migMatNameProp
} else {
migText="nomig"
curCex=migMatNameProp/2
slide=slide*2
}
text(0+slide, time2DrawArrows, labels=migText, cex=curCex, col=migrMatCol)
lastSink=he.sink
}
}
#-- Draw last branch
segments(activePops[1],topOfCircle[activePops[1]], activePops[1], yTimeLimit)
#============================== PLOT LEGENDS IN MARGINS ========================
#Compute space available in margin
minY.coo=grconvertY(0, from="nic", to="user")
par(xpd=NA)
#--- Draw population size scale with circles of different sizes .................
maxOrder=ceiling(log10(maxPopSize))
minOrder=floor(log10(minPopSize))
popSizeRadius=10^(maxOrder:minOrder)
winWidth=numSamples+2
ypos=3/4*minY.coo
text(x=-winWidth/10*1.2, y=ypos, labels="Pop. \nsizes ", cex=.8, pos=2)
for (i in 1:length(popSizeRadius)) {
curRadius=interpolRadius(popSizeRadius[i], minPopSize, maxPopSize, minRadius, maxRadius, drawLogPopSize)
# print(curRadius)
if (curRadius>0) {
xpos=-winWidth/10+(i-1)*winWidth/10
draw.circle(xpos, ypos, radius=curRadius, col=popCol, border=popBorderCol)
}
text(xpos, ypos-abs(ypos)*0.1, popSizeRadius[i], cex=0.7, pos=1, col="black")
}
#--- Legend for growing or shrinking populations ...............................
if (isGrowth) {
x=winWidth-1.5*winWidth/10; y=ypos+abs(ypos)*0.2
text(x, y-abs(ypos)*0.1, labels="Populations", cex=0.8)
x=winWidth-2*winWidth/10; y=ypos-abs(ypos)*0.1
drawTriangle(1, x, y, size=0.15, aspRatio, color=growthCol)
text(x, y-abs(ypos)*0.1, labels="growing", pos=NULL, cex=0.7)
x=winWidth-winWidth/10
drawTriangle(-1, x, y, size=0.15, aspRatio, color=growthCol)
text(x, y-abs(ypos)*0.1, labels="shrinking", pos=NULL, cex=0.7)
}
#if (printPDF) dev.off()
dev.off()
## null device
## 1
The same process described for V. destructor demographic history inferences was independently applied for its sister species V. jacobsoni.
Contrary to V. destructor for which the scenario 6 was undeniably the most likely according to both AIC and model likelihood distribution, both scenarios 3 and 6 seems likely in the case of recent host switch by V. jacobsoni.
As a consequence, we repeated the same process of bootstraps computations as done for V. destructor but for both scenarios: DIV + BOT and IM + BOT.
vj.lhoods1050 <- read.csv("fsc_run/lhoods_1050SNPS_vjpng.csv", header = TRUE)
vj.lhoods10500 <- read.csv("fsc_run/lhoods_10500SNPS_vjpng.csv", header = TRUE)
vj.lhoods35000 <- read.csv("fsc_run/lhoods_35000SNPS_vjpng.csv", header = TRUE)
vj.lhoods105000 <- read.csv("fsc_run/lhoods_105000SNPS_vjpng.csv", header = TRUE)
par(mfrow=c(1,4))
#Plot 1 for likelihoods with smaller subset
boxplot(range = 0,vj.lhoods1050$scenario1_div [3:102], vj.lhoods1050$scenario2_divGwt2n[3:102], vj.lhoods1050$scenario3_divGwt[3:102], vj.lhoods1050$scenario4_IM[3:102], vj.lhoods1050$scenario5_IMGwt2n[3:102], vj.lhoods1050$scenario6_IMGwt[3:102], ylab="log 10 Likelihood",xaxt="n")
points(1,vj.lhoods1050$scenario1_div [2],col = "blue", pch = 20)
points(2,vj.lhoods1050$scenario2_divGwt2n [2],col = "blue", pch = 20)
points(3,vj.lhoods1050$scenario3_divGwt [2],col = "blue", pch = 20)
points(4,vj.lhoods1050$scenario4_IM [2],col = "blue", pch = 20)
points(5,vj.lhoods1050$scenario5_IMGwt2n [2],col = "blue", pch = 20)
points(6,vj.lhoods1050$scenario6_IMGwt [2],col = "blue", pch = 20)
axis(side=1,at=1:6, labels=c("DIV1","DIV2","DIV3","IM4","IM5","IM6"))
title("V. jacobsoni 472 SNPS")
#Plot 2 for likelihoods with second smaller subset
boxplot(range = 0,vj.lhoods10500$scenario1_div [3:102], vj.lhoods10500$scenario2_divGwt2n[3:102], vj.lhoods10500$scenario3_divGwt[3:102], vj.lhoods10500$scenario4_IM[3:102], vj.lhoods10500$scenario5_IMGwt2n[3:102], vj.lhoods10500$scenario6_IMGwt[3:102], ylab="log 10 Likelihood",xaxt="n")
points(1,vj.lhoods10500$scenario1_div [2],col = "blue", pch = 20)
points(2,vj.lhoods10500$scenario2_divGwt2n [2],col = "blue", pch = 20)
points(3,vj.lhoods10500$scenario3_divGwt [2],col = "blue", pch = 20)
points(4,vj.lhoods10500$scenario4_IM [2],col = "blue", pch = 20)
points(5,vj.lhoods10500$scenario5_IMGwt2n [2],col = "blue", pch = 20)
points(6,vj.lhoods10500$scenario6_IMGwt [2],col = "blue", pch = 20)
axis(side=1,at=1:6, labels=c("DIV1","DIV2","DIV3","IM4","IM5","IM6"))
title("V. jacobsoni 4,713 SNPS")
#Plot 3 for likelihoods with third subset
boxplot(range = 0,vj.lhoods35000$scenario1_div [3:102], vj.lhoods35000$scenario2_divGwt2n[3:102], vj.lhoods35000$scenario3_divGwt[3:102], vj.lhoods35000$scenario4_IM[3:102], vj.lhoods35000$scenario5_IMGwt2n[3:102], vj.lhoods35000$scenario6_IMGwt[3:102], ylab="log 10 Likelihood",xaxt="n")
points(1,vj.lhoods35000$scenario1_div [2],col = "blue", pch = 20)
points(2,vj.lhoods35000$scenario2_divGwt2n [2],col = "blue", pch = 20)
points(3,vj.lhoods35000$scenario3_divGwt [2],col = "blue", pch = 20)
points(4,vj.lhoods35000$scenario4_IM [2],col = "blue", pch = 20)
points(5,vj.lhoods35000$scenario5_IMGwt2n [2],col = "blue", pch = 20)
points(6,vj.lhoods35000$scenario6_IMGwt [2],col = "blue", pch = 20)
axis(side=1,at=1:6, labels=c("DIV1","DIV2","DIV3","IM4","IM5","IM6"))
title("V. jacobsoni 15,363 SNPS")
#Plot 4 for likelihoods with larger subset
boxplot(range = 0,vj.lhoods105000$scenario1_div [3:102], vj.lhoods105000$scenario2_divGwt2n[3:102], vj.lhoods105000$scenario3_divGwt[3:102], vj.lhoods105000$scenario4_IM[3:102], vj.lhoods105000$scenario5_IMGwt2n[3:102], vj.lhoods105000$scenario6_IMGwt[3:102], ylab="log 10 Likelihood",xaxt="n")
points(1,vj.lhoods105000$scenario1_div [2],col = "blue", pch = 20)
points(2,vj.lhoods105000$scenario2_divGwt2n [2],col = "blue", pch = 20)
points(3,vj.lhoods105000$scenario3_divGwt [2],col = "blue", pch = 20)
points(4,vj.lhoods105000$scenario4_IM [2],col = "blue", pch = 20)
points(5,vj.lhoods105000$scenario5_IMGwt2n [2],col = "blue", pch = 20)
points(6,vj.lhoods105000$scenario6_IMGwt [2],col = "blue", pch = 20)
axis(side=1,at=1:6, labels=c("DIV1","DIV2","DIV3","IM4","IM5","IM6"))
title("V. jacobsoni 46,416 SNPS")
# import the file containing the 100 best runs
vj.boot3 <- read.csv("fsc_run/png_105000nopv_scen3_BOOT.csv", header = TRUE)
# Effective size for Varroa mites on original host Apis cerana
N_VjAc1.boot = boot(vj.boot3$NVJAC1, my.mean, 10000)
my.mean(vj.boot3$NVJAC1,1:length(vj.boot3$NVJAC1))
## [1] 228.34
boot.ci(N_VjAc1.boot)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = N_VjAc1.boot)
##
## Intervals :
## Level Normal Basic
## 95% (208.7, 248.0 ) (208.0, 247.6 )
##
## Level Percentile BCa
## 95% (209.1, 248.6 ) (209.5, 249.0 )
## Calculations and Intervals on Original Scale
# Effective size for Varroa mites on new host Apis mellifera
N_VjAm0.boot = boot(vj.boot3$NVJAM0, my.mean, 10000)
my.mean(vj.boot3$NVJAM0,1:length(vj.boot3$NVJAM0))
## [1] 6708.25
boot.ci(N_VjAm0.boot)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = N_VjAm0.boot)
##
## Intervals :
## Level Normal Basic
## 95% (6231, 7193 ) (6229, 7191 )
##
## Level Percentile BCa
## 95% (6226, 7188 ) (6235, 7198 )
## Calculations and Intervals on Original Scale
# Number of generations when new population of Varroa mites was founded
TBOT.boot = boot(vj.boot3$TBOT, my.mean, 10000)
my.mean(vj.boot3$TBOT,1:length(vj.boot3$TBOT))
## [1] 58.66
boot.ci(TBOT.boot)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = TBOT.boot)
##
## Intervals :
## Level Normal Basic
## 95% (53.44, 63.74 ) (53.34, 63.67 )
##
## Level Percentile BCa
## 95% (53.65, 63.98 ) (53.67, 64.01 )
## Calculations and Intervals on Original Scale
# Growth rate of Varroa mites on new host Apis mellifera
GAM.boot = boot(vj.boot3$GAM, my.mean, 10000)
my.mean(vj.boot3$GAM,1:length(vj.boot3$GAM))
## [1] -0.156716
boot.ci(GAM.boot)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = GAM.boot)
##
## Intervals :
## Level Normal Basic
## 95% (-0.1730, -0.1405 ) (-0.1722, -0.1400 )
##
## Level Percentile BCa
## 95% (-0.1734, -0.1413 ) (-0.1750, -0.1423 )
## Calculations and Intervals on Original Scale
# Effective size for Varroa mites able to jump on A. mellifera to create new population
JUMP.boot = boot(vj.boot3$JUMP, my.mean, 10000)
my.mean(vj.boot3$JUMP,1:length(vj.boot3$JUMP))
## [1] 348.86
boot.ci(JUMP.boot)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = JUMP.boot)
##
## Intervals :
## Level Normal Basic
## 95% (313.9, 384.3 ) (313.2, 382.8 )
##
## Level Percentile BCa
## 95% (314.9, 384.5 ) (317.2, 388.7 )
## Calculations and Intervals on Original Scale
#...........................................................................................
# (c) Laurent Excoffier and Vitor Sousa June-November 2015
#
# Small R program to draw the evolutonary scenario described by a given par file
# This is mainly for visual checking that the modeled scenario corresponds to
# what was intended
#
#...........................................................................................
#
args=commandArgs(TRUE)
print(args)
## character(0)
#Expects par file name on command line
if(length(args)) {
parFileName=args[1]
} else {
#REPLACE BY THE NAME OF THE PAR FILE YOU WANT TO ANALYSE
parFileName="fsc_run/VJM3_downsample_divGwt_29_boot.par"
}
separator=" "
migrMatCol="coral"
admixCol="blue"
popFusionColor="black"
popCol="lightgrey"
popBorderCol="black"
inbreedColor="lightblue"
oldSampColor="olivedrab4"
timeCol="tan4"
growthCol="hotpink3"
propLastsegment=0.05
migMatNameProp=0.8
migMatLineLength=0.3
timeProp=0.6
maxRadius=1/40
minRadius=maxRadius/3
arrowLength=0.2
logScaleAxis=""
timeOffset=0.25
migrOffset=0.1
curvedArrowLTY=1
drawLogPopSize=T
plotMigrRates=T
migrRateTextSizeCEX=0.5
#Define plot area size for PDF
pdf.x.size=7
pdf.y.size=10
rescalingFactor=1.0 #Don't touch this!
printPDF=F
# parFile=readLines(con=parFileName)
########################### Reading par file #########################
parFile=scan(parFileName, character(0), sep = "\n", strip.white = TRUE) # separate each line
pdfFileName=paste(parFileName, ".pdf", sep="")
#--- Clean input file by removing consecutive separators, and keep ---------------
#--- Function to remove separators within a string
removeTrailingSep=function(string, sep) {
temp=strsplit(string, split=sep)
temp2=temp[[1]][nchar(temp[[1]])>0]
cleanStr=temp2[1]
if (length(temp2)>1) {
for (i in 2:length(temp2)) {
cleanStr=paste(cleanStr, temp2[i], sep=sep)
}
}
cleanStr
}
#--- Replace Keep by -9999
replaceKeep=function(string) {
if (grepl("keep", string)) {
return(gsub("keep", "-9999", string))
}
return(string)
}
#Remove both multiple consecutive whitespace and tabs and replace keep keyword
for (i in 1:length(parFile)) {
parFile[i]=removeTrailingSep(parFile[i], sep='\t')
parFile[i]=removeTrailingSep(parFile[i], sep=' ')
parFile[i]=replaceKeep(parFile[i])
}
#-------------------------------------------------------------------------------
#--- Get number of samples on line 2 -----
l.numsamples=parFile[2]
sp.l=unlist(strsplit(l.numsamples, split=' '))
tab.l=unlist(strsplit(l.numsamples, split='\t'))
if (length(sp.l)>=length(tab.l)) {
numSamples=as.numeric(sp.l[1])
separator=" "
} else {
numSamples=as.numeric(tab.l[1])
separator="\t"
}
#--- Reading numbers on separate lines -----
getNumbers=function(start, parFile, numSamples) {
for (i in 1:numSamples) {
curnum=as.numeric(unlist(strsplit(parFile[start+i], split=separator))[1])
if (i==1) {
num=curnum
} else {
num=c(num, curnum)
}
}
num
}
#--- Get population sizes -----------------
start=3
popSizes=getNumbers(start, parFile, numSamples)
#Rescaling pop sizes
popSizes=round(popSizes*rescalingFactor, digits=0)
iniPopSizes=popSizes
#--- Get sample sizes -----------------
readSampleSizesTimesAndInbreedingLevel=function(start, parFile, numSamples) {
for (i in 1:numSamples) {
curLine=unlist(strsplit(parFile[start+i], split=separator))
curSampSize=as.numeric(curLine[1])
curSampTime=0
curInbreeding=0
if (length(curLine)>1) curSampTime=as.numeric(curLine[2])
if (length(curLine)>2) curInbreeding=as.numeric(curLine[3])
if (i==1) {
sampSize=curSampSize
sampTime=curSampTime
inbreeding=curInbreeding
} else {
sampSize=c(sampSize,curSampSize)
sampTime=c(sampTime,curSampTime)
inbreeding=c(inbreeding,curInbreeding)
}
}
list(ss=sampSize, st=sampTime, inb=inbreeding)
}
start=start+numSamples+1
# sampSizes=getNumbers(start, parFile, numSamples)
sampSizesStats=readSampleSizesTimesAndInbreedingLevel(start, parFile, numSamples)
#Rescaling sample times
sampSizesStats$st=round(sampSizesStats$st*rescalingFactor, digits=0)
sampSizes=sampSizesStats$ss
sampTimes=sampSizesStats$st
inbrCoeff=sampSizesStats$inb
#--- Get growth rates -----------------
start=start+numSamples+1
growthRatesInitial=getNumbers(start, parFile, numSamples)
# save this into growthRates which will be used and updated when printing historical events
growthRates=growthRatesInitial
#--- Get number of migration matrices -----------------
start=start+numSamples+1
numMigMat=as.numeric(unlist(strsplit(parFile[start+1], split=separator))[1])
#--- Read migration matrix
readMigMat=function(start, parFile, numSamples) {
for (i in 1:numSamples) {
curmigs=as.numeric(unlist(strsplit(parFile[start+i], split=separator)))
if (i==1) {
migs=curmigs
} else {
migs=rbind(migs, curmigs)
}
}
rownames(migs)=1:numSamples
migs
}
#--- Get migration matrices as a list --------------
start=start+2
migMats=list()
if (numMigMat>0) {
for (i in 1:numMigMat) {
curMigMat=readMigMat(start, parFile, numSamples)
migMats[[i]]=curMigMat
start=start+numSamples+1
}
}
#Rescaling migration rates
if (numMigMat>0) {
for (i in 1:numMigMat) {
migMats[[i]]=migMats[[i]]/rescalingFactor
}
}
#--- Get number of historical events
start=start+1
numHistEvents=as.numeric(unlist(strsplit(parFile[start], split=separator))[1])
###################### HISTORICAL EVENTS HANDLING ##############################
#..... Read Historical Event .......
last.he.time=0
if (numHistEvents>0) {
for (i in 1:numHistEvents) {
start=start+1
#Take care of nomig keyword
nomig=F
if (grepl("nomig", parFile[start])) {
nomig=T
gsub("nomig", "", parFile[start])
}
curHE=as.numeric(unlist(strsplit(parFile[start], split=separator)))
if (nomig) curHE[7]=-1
#Rescaling time of event
curHE[1]=round(curHE[1]*rescalingFactor, digits=0)
if (i==1) {
histEvents=curHE
last.he.time=curHE[1]
} else {
histEvents=rbind(histEvents, curHE)
if (histEvents[i,1] > last.he.time) last.he.time=histEvents[i,1]
}
}
if (numHistEvents>1) rownames(histEvents)=1:numHistEvents
}
names(last.he.time)=""
last.he.time=as.numeric(last.he.time)
yTimeLimit=0
if (last.he.time!=0) yTimeLimit=last.he.time*(1+propLastsegment) #Add propLastsegment to y axis after last event (to draw stuff)
#Reorder events by their times
if (numHistEvents>1) histEvents=histEvents[order(histEvents[,1],decreasing=FALSE),] else {
if (numHistEvents==1) histEvents=matrix(histEvents, nrow=1, byrow = T)
}
endReadParFile=start
############################## PLOTTING THE MODELED SCENARIO ######################################
#--- Graphical functions ...........................................................................
fullHeadArrow=function(x0, y0, x1, y1, length, angle, color="black", weight=1) {
arrows(x0, y0, x1, y1, length, angle, code=2, lty=1, col=color, lwd=weight)
arrows(x0, y0, x1, y1, length, angle*0.80, code=2, lty=1, col=color, lwd=weight)
arrows(x0, y0, x1, y1, length, angle*0.60, code=2, lty=1, col=color, lwd=weight)
arrows(x0, y0, x1, y1, length, angle*0.40, code=2, lty=1, col=color, lwd=weight)
arrows(x0, y0, x1, y1, length, angle*0.20, code=2, lty=1, col=color, lwd=weight)
arrows(x0, y0, x1, y1, length, angle*0.10, code=2, lty=1, col=color, lwd=weight)
}
drawTriangle=function(growth, x, y, size, aspRatio, color) {
if (growth>0) {
x0=x; y0=y
x1=x-size/2; y1=y+size/2*aspRatio
x2=x+size/2; y2=y+size/2*aspRatio
} else {
x0=x-size/2; y0=y
x1=x; y1=y+size/2*aspRatio
x2=x+size/2; y2=y
}
polygon(c(x0, x1, x2, x0), c(y0, y1, y2, y0), col=color)
return(y+size/2*aspRatio)
}
#--- Computing maximum current pop size ............................................................
maxPopSize=popSizes[1]
minPopSize=popSizes[1]
if (length(popSizes)>1){
for (i in 2:length(popSizes)) {
if (popSizes[i]>maxPopSize) maxPopSize=popSizes[i];
if (popSizes[i]<minPopSize) minPopSize=popSizes[i];
}
}
#--- Find min and max pop sizes over the whole population history ..................................
ps=popSizes
isGrowth=FALSE
if (numHistEvents) {
#Need to keep track ofgrowth rates over time
gRates=growthRates
prevTime=0
for (i in 1:numHistEvents) {
he=histEvents[i,]
curTime=he[1]
sink=he[3]+1
resize=he[5]
growth=he[6]
#Begin by resizing pop sizes due to growth since last event
for (j in 1:length(ps)) {
ps[j]=ps[j]*exp(gRates[j]*(curTime-prevTime))
if (ps[j]>maxPopSize) maxPopSize=ps[j]
if (ps[j]<minPopSize) minPopSize=ps[j]
if (ps[j]==0) print(paste("Deme ", j-1, " reaches size zero at time ", curTime, " due to large negative growth (", gRates[j], ")", sep=""))
if (is.infinite(ps[j])) print(paste("Deme ", j-1, " reaches infinite size at time ", curTime, " due to large positive growth (", gRates[j], ")", sep=""))
}
prevTime=curTime
#Handle transformed keep keyword
if(growth!=-9999) {
gRates[sink]=growth
}
#Implement resize
ps[sink]=ps[sink]*resize
if (ps[sink]>maxPopSize) maxPopSize=ps[sink]
if (ps[sink]<minPopSize) minPopSize=ps[sink]
}
}
#-- Function to compute the circle radius for a given pop size ....................................
interpolRadius=function(curSize, minSize, maxSize, minRadius, maxRadius, logScale) {
if(logScale) {
minSize=log10(minSize)
maxSize=log10(maxSize)
curSize=log10(curSize)
}
curRadius=minRadius+(curSize-minSize)*(maxRadius-minRadius)/(maxSize-minSize)
curRadius
}
#============================== ========================
#============================== BEGIN MAIN PLOT ========================
#============================== ========================
library("plotrix")
library("diagram")
if (printPDF) {
pdf(pdfFileName, width=pdf.x.size, height=pdf.y.size)
}
par(xpd=F, mar=c(6,4,3,0.5))
maxRadius=maxRadius*(numSamples+2)
minRadius=minRadius*(numSamples+2)
title=parFileName
plot(x=1:numSamples, type="n", xlab="", ylab="time (gen)", xlim=c(-0.5, numSamples+1.5), cex.main=0.8,
ylim=c(0, yTimeLimit), main=title, xaxt = 'n', log=logScaleAxis, cex.axis=0.9, cex.lab=0.9)
axis(side=1, labels=c("Mig Mat",0:(numSamples-1), " \nTimes"), at=0:(numSamples+1), cex.axis=0.8)
mtext("demes", side=1,line=2, cex=0.8)
w <- par("pin")[1]/diff(par("usr")[1:2])
h <- par("pin")[2]/diff(par("usr")[3:4])
aspRatio <- w/h
# aspRatio=yTimeLimit/(numSamples+2)
#--- Plotting population circles according to their current size and sampling time
for (i in 1:numSamples) {
curRadius=interpolRadius(popSizes[i], minPopSize, maxPopSize, minRadius, maxRadius, drawLogPopSize)
# print(curRadius)
if (i==1) topOfCircle=sampTimes[i]+curRadius*aspRatio else topOfCircle=c(topOfCircle, sampTimes[i]+curRadius*aspRatio)
curColor=popCol
if (sampSizes[i]==0) curColor="white"
if (inbrCoeff[i]==0) {
draw.circle(i, sampTimes[i], radius = curRadius, col=curColor, border=popBorderCol)
} else {
floating.pie(i, sampTimes[i], c(inbrCoeff[i], 1-inbrCoeff[i]), radius = curRadius,
col=c(inbreedColor,curColor), startpos=pi, border=popBorderCol)
}
if (sampTimes[i]>0) {
text(i, sampTimes[i]-curRadius*aspRatio, labels=sampTimes[i], cex=timeProp, col=oldSampColor, pos=1)
}
#Draw a vertical arrow in case of pop growth
curGrowthRate=growthRates[i]
#--- Draw gtrowing or shrinking triangle on top of pop circle to reflect growth type
if (curGrowthRate!=0) {
arLength=0.15
topOfCircle[i]=drawTriangle(curGrowthRate, x=i, y=topOfCircle[i], arLength, aspRatio, color=growthCol)
}
}
#--- Handle first migration matrix .......................
curMigMatNum=0
curvature=0.0075*last.he.time
text(0-migrOffset, 0, labels=0, cex=migMatNameProp, col=migrMatCol)
if (numMigMat) {
curMigMat=migMats[1][[1]]
for (sink in 1:numSamples) {
for (sourc in 1:numSamples) {
if (sink!=sourc & curMigMat[sourc, sink]>0) {
differ=sourc-sink
curvedarrow(from=c(sourc, 0), to=c(sink,0), curve=curvature*(abs(differ)*0.55^abs(differ)), arr.adj=1,
arr.pos=0.5, arr.type="triangle", arr.col=migrMatCol, lwd=1, lty=curvedArrowLTY,
lcol=migrMatCol, arr.length=arrowLength)
if (plotMigrRates) {
curNm=round(curMigMat[sourc, sink]*popSizes[sourc], digits=2)
if (differ>0) {
text(sink+abs(differ)/2, aspRatio*0.15*abs(differ), labels=curNm, cex=migrRateTextSizeCEX, col=migrMatCol)
} else {
text(sourc+abs(differ)/2, -aspRatio*0.15*abs(differ), labels=curNm, cex=migrRateTextSizeCEX, col=migrMatCol)
}
}
}
}
}
}
#---Draw all events on the population tree
lastTime=0
activePops=1:numSamples
numActivePops=numSamples
lastSink=-1
if (numHistEvents) {
for (i in 1:numHistEvents) {
#Extract historical event
he=histEvents[i,]
he.time=he[1]
he.source=he[2]+1 #+1 is due to the use of base 0 for deme number in fsc
he.sink=he[3]+1
he.migr=he[4]
he.resize=he[5]
he.growth=he[6]
if(he.growth==-9999) { #handle transformed keep keyword
he.growth=growthRates[he.sink] #Keep current growth rate
}
he.migrMat=he[7]
if(he.migrMat==-9999) { #handle transformed keep keyword
he.migrMat=curMigMatNum
}
#-- Draw event time
if (i%%2) {
slide=timeOffset
}
else {
slide=-timeOffset
}
#Draw all vertical segments ....................
for (p in 1:length(activePops)) {
#Handle population with older sampling times
if (sampTimes[activePops[p]]>0) {
minTime=max(lastTime, topOfCircle[activePops[p]]);
} else {
minTime=lastTime;
}
if (he.time>sampTimes[activePops[p]] & he.time>topOfCircle[activePops[p]]) {
if (activePops[p]==lastSink | i==1) {
segments(activePops[p], topOfCircle[activePops[p]], activePops[p], he.time)
} else {
segments(activePops[p], minTime, activePops[p], he.time)
}
}
}
#Handle growth rate changes since last event ......................
#Update pop sizes according to current growth rates
if (i>1) {
prev.he=histEvents[i-1,]
branchLength=he.time-lastTime
for (p in 1:length(activePops)) {
curPop=activePops[p]
popSizes[curPop]=popSizes[curPop]*exp(growthRates[curPop]*branchLength);
}
}
#Update growth rate
growthRates[he.sink]=he.growth
lastTime=he.time
#Handle resize of sink pop ........................
popSizes[he.sink]=popSizes[he.sink]*he.resize
curRadius=interpolRadius(popSizes[he.sink], minPopSize, maxPopSize, minRadius, maxRadius, drawLogPopSize)
# popRadius[he.sink]=curRadius
topOfCircle[he.sink]=he.time+curRadius*aspRatio
draw.circle(he.sink, he.time, radius = curRadius, col=popCol, border=popBorderCol)
#--- Draw growing or shrinking triangle on top of pop circle to reflect growth type
if (he.growth!=0) {
arLength=0.15
topOfCircle[he.sink]=drawTriangle(he.growth, x=he.sink, y=topOfCircle[he.sink], arLength, aspRatio, color=growthCol)
}
#Handle population fusion .........................
if (he.migr>=1 & he.sink!=he.source) { #This is a population fusion
if (numActivePops==numSamples) removedPops=he.source else removedPops=c(removedPops, he.source)
numActivePops=numActivePops-1
activePops=(1:numSamples)[-removedPops]
#Draw connecting arrows from source to sink
fullHeadArrow(he.source, he.time, he.sink, he.time, length=0.15, angle=20)
#Redraw time with the right color
text(numSamples+1+slide, he.time, labels=he.time, cex=timeProp, col=popFusionColor)
} else {
#--- Handle admixture event ........................
if (he.migr>0 & he.migr<1) {
#Draw connecting arrows from source to sink
segments(he.source, he.time, he.sink, he.time, col=admixCol, lty=2)
if (he.sink>he.source) {
fullHeadArrow(he.sink-0.15, he.time, he.sink, he.time, length=0.15, angle=20, color=admixCol)
} else {
fullHeadArrow(he.sink+0.15, he.time, he.sink, he.time, length=0.15, angle=20, color=admixCol)
}
#Redraw time with the right color
text(numSamples+1+slide, he.time, labels=he.time, cex=timeProp, col=admixCol)
}
else text(numSamples+1+slide, he.time, labels=he.time, cex=timeProp, col=timeCol)
}
#--- Handle migmat change ............................
if (i!=numHistEvents) nextTime=histEvents[i+1,1] else {
nextTime=yTimeLimit
}
time2DrawArrows=(he.time+nextTime)/2
if (he.migrMat!=curMigMatNum) {
if (he.migrMat>-1) {
curMigMat=migMats[he.migrMat+1][[1]]
for (sink in 1:numSamples) {
for (sourc in 1:numSamples) {
if (sink!=sourc & curMigMat[sourc, sink]>0) {
differ=sourc-sink
curvedarrow(from=c(sourc, time2DrawArrows), to=c(sink, time2DrawArrows), curve=curvature*(abs(differ)*0.55^abs(differ)),
arr.adj=1, arr.pos=0.5, arr.type="triangle", arr.col=migrMatCol, lwd=1, lty=curvedArrowLTY,
lcol=migrMatCol, arr.length=arrowLength)
if (plotMigrRates) {
#Write Nm values
curNm=round(curMigMat[sourc, sink]*popSizes[sourc], digits=2)
if (differ>0) {
text(sink+abs(differ)/2, time2DrawArrows+aspRatio*0.15*abs(differ), labels=curNm, cex=migrRateTextSizeCEX, col=migrMatCol)
} else {
text(sourc+abs(differ)/2, time2DrawArrows-aspRatio*0.15*abs(differ), labels=curNm, cex=migrRateTextSizeCEX, col=migrMatCol)
}
}
}
}
# curMigMatNum=he.migrMat
}
}
#--- Draw separation between migration matrices numbers on the left
segments(-migMatLineLength/2, he.time, migMatLineLength/2, he.time, lty=3, col=migrMatCol)
}
curMigMatNum=he.migrMat
#Output current valid migration matrix
if (i%%2) {
slide=migrOffset
}
else {
slide=-migrOffset
}
if (he.migrMat>-1) {
migText=he.migrMat
curCex=migMatNameProp
} else {
migText="nomig"
curCex=migMatNameProp/2
slide=slide*2
}
text(0+slide, time2DrawArrows, labels=migText, cex=curCex, col=migrMatCol)
lastSink=he.sink
}
}
#-- Draw last branch
segments(activePops[1],topOfCircle[activePops[1]], activePops[1], yTimeLimit)
#============================== PLOT LEGENDS IN MARGINS ========================
#Compute space available in margin
minY.coo=grconvertY(0, from="nic", to="user")
par(xpd=NA)
#--- Draw population size scale with circles of different sizes .................
maxOrder=ceiling(log10(maxPopSize))
minOrder=floor(log10(minPopSize))
popSizeRadius=10^(maxOrder:minOrder)
winWidth=numSamples+2
ypos=3/4*minY.coo
text(x=-winWidth/10*1.2, y=ypos, labels="Pop. \nsizes ", cex=.8, pos=2)
for (i in 1:length(popSizeRadius)) {
curRadius=interpolRadius(popSizeRadius[i], minPopSize, maxPopSize, minRadius, maxRadius, drawLogPopSize)
# print(curRadius)
if (curRadius>0) {
xpos=-winWidth/10+(i-1)*winWidth/10
draw.circle(xpos, ypos, radius=curRadius, col=popCol, border=popBorderCol)
}
text(xpos, ypos-abs(ypos)*0.1, popSizeRadius[i], cex=0.7, pos=1, col="black")
}
#--- Legend for growing or shrinking populations ...............................
if (isGrowth) {
x=winWidth-1.5*winWidth/10; y=ypos+abs(ypos)*0.2
text(x, y-abs(ypos)*0.1, labels="Populations", cex=0.8)
x=winWidth-2*winWidth/10; y=ypos-abs(ypos)*0.1
drawTriangle(1, x, y, size=0.15, aspRatio, color=growthCol)
text(x, y-abs(ypos)*0.1, labels="growing", pos=NULL, cex=0.7)
x=winWidth-winWidth/10
drawTriangle(-1, x, y, size=0.15, aspRatio, color=growthCol)
text(x, y-abs(ypos)*0.1, labels="shrinking", pos=NULL, cex=0.7)
}
#if (printPDF) dev.off()
dev.off()
## null device
## 1
# import the file containing the 100 best runs
vj.boot6 <- read.csv("fsc_run/png_105000nopv_scen6_BOOT.csv", header = TRUE)
# Effective size for Varroa mites on original host Apis cerana
N_VjAc1.boot = boot(vj.boot6$NVJAC1, my.mean, 10000)
my.mean(vj.boot6$NVJAC1,1:length(vj.boot6$NVJAC1))
## [1] 150.5
boot.ci(N_VjAc1.boot)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = N_VjAc1.boot)
##
## Intervals :
## Level Normal Basic
## 95% (137.2, 163.8 ) (137.0, 163.4 )
##
## Level Percentile BCa
## 95% (137.6, 164.0 ) (138.5, 165.0 )
## Calculations and Intervals on Original Scale
# Effective size for Varroa mites on new host Apis mellifera
N_VjAm0.boot = boot(vj.boot6$NVJAM0, my.mean, 10000)
my.mean(vj.boot6$NVJAM0,1:length(vj.boot6$NVJAM0))
## [1] 7922.38
boot.ci(N_VjAm0.boot)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = N_VjAm0.boot)
##
## Intervals :
## Level Normal Basic
## 95% (7406, 8435 ) (7398, 8431 )
##
## Level Percentile BCa
## 95% (7414, 8447 ) (7427, 8459 )
## Calculations and Intervals on Original Scale
# Number of generations when new population of Varroa mites was founded
TBOT.boot = boot(vj.boot6$TBOT, my.mean, 10000)
my.mean(vj.boot6$TBOT,1:length(vj.boot6$TBOT))
## [1] 38.99
boot.ci(TBOT.boot)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = TBOT.boot)
##
## Intervals :
## Level Normal Basic
## 95% (35.52, 42.43 ) (35.41, 42.28 )
##
## Level Percentile BCa
## 95% (35.70, 42.57 ) (35.94, 42.85 )
## Calculations and Intervals on Original Scale
# Growth rate of Varroa mites on new host Apis mellifera
GAM.boot = boot(vj.boot6$GAM, my.mean, 10000)
my.mean(vj.boot6$GAM,1:length(vj.boot6$GAM))
## [1] -0.253887
boot.ci(GAM.boot)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = GAM.boot)
##
## Intervals :
## Level Normal Basic
## 95% (-0.2776, -0.2302 ) (-0.2770, -0.2292 )
##
## Level Percentile BCa
## 95% (-0.2786, -0.2308 ) (-0.2799, -0.2320 )
## Calculations and Intervals on Original Scale
# Effective size for Varroa mites able to jump on A. mellifera to create new population
JUMP.boot = boot(vj.boot6$JUMP, my.mean, 10000)
my.mean(vj.boot6$JUMP,1:length(vj.boot6$JUMP))
## [1] 295.06
boot.ci(JUMP.boot)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = JUMP.boot)
##
## Intervals :
## Level Normal Basic
## 95% (268.8, 321.3 ) (268.8, 320.3 )
##
## Level Percentile BCa
## 95% (269.8, 321.4 ) (271.3, 323.0 )
## Calculations and Intervals on Original Scale
N0M1.boot = boot(vj.boot6$N0M1, my.mean, 10000)
my.mean(vj.boot6$N0M1,1:length(vj.boot6$N0M1))
## [1] 0.02610133
boot.ci(N0M1.boot)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = N0M1.boot)
##
## Intervals :
## Level Normal Basic
## 95% ( 0.0166, 0.0354 ) ( 0.0160, 0.0348 )
##
## Level Percentile BCa
## 95% ( 0.0174, 0.0362 ) ( 0.0183, 0.0376 )
## Calculations and Intervals on Original Scale
N1M0.boot = boot(vj.boot6$N1M0, my.mean, 10000)
my.mean(vj.boot6$N1M0,1:length(vj.boot6$N1M0))
## [1] 0.00691537
boot.ci(N1M0.boot)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = N1M0.boot)
##
## Intervals :
## Level Normal Basic
## 95% ( 0.0044, 0.0094 ) ( 0.0044, 0.0093 )
##
## Level Percentile BCa
## 95% ( 0.0045, 0.0095 ) ( 0.0048, 0.0098 )
## Calculations and Intervals on Original Scale
MACtoAM.boot = boot(vj.boot6$MACtoAM, my.mean, 10000)
my.mean(vj.boot6$MACtoAM,1:length(vj.boot6$MACtoAM))
## [1] 5.264404e-05
boot.ci(MACtoAM.boot)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = MACtoAM.boot)
##
## Intervals :
## Level Normal Basic
## 95% ( 0.0000, 0.0001 ) ( 0.0000, 0.0001 )
##
## Level Percentile BCa
## 95% ( 0.0000, 0.0001 ) ( 0.0000, 0.0001 )
## Calculations and Intervals on Original Scale
MAMtoAC.boot = boot(vj.boot6$MAMtoAC, my.mean, 10000)
my.mean(vj.boot6$MAMtoAC,1:length(vj.boot6$MAMtoAC))
## [1] 3.453322e-06
boot.ci(MAMtoAC.boot)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = MAMtoAC.boot)
##
## Intervals :
## Level Normal Basic
## 95% ( 0, 0 ) ( 0, 0 )
##
## Level Percentile BCa
## 95% ( 0, 0 ) ( 0, 0 )
## Calculations and Intervals on Original Scale
#...........................................................................................
# (c) Laurent Excoffier and Vitor Sousa June-November 2015
#
# Small R program to draw the evolutonary scenario described by a given par file
# This is mainly for visual checking that the modeled scenario corresponds to
# what was intended
#
#...........................................................................................
#
args=commandArgs(TRUE)
print(args)
## character(0)
#Expects par file name on command line
if(length(args)) {
parFileName=args[1]
} else {
#REPLACE BY THE NAME OF THE PAR FILE YOU WANT TO ANALYSE
parFileName="fsc_run/VJM6_downsample_IMGwt_32_boot.par"
}
separator=" "
migrMatCol="coral"
admixCol="blue"
popFusionColor="black"
popCol="lightgrey"
popBorderCol="black"
inbreedColor="lightblue"
oldSampColor="olivedrab4"
timeCol="tan4"
growthCol="hotpink3"
propLastsegment=0.05
migMatNameProp=0.8
migMatLineLength=0.3
timeProp=0.6
maxRadius=1/40
minRadius=maxRadius/3
arrowLength=0.2
logScaleAxis=""
timeOffset=0.25
migrOffset=0.1
curvedArrowLTY=1
drawLogPopSize=T
plotMigrRates=T
migrRateTextSizeCEX=0.5
#Define plot area size for PDF
pdf.x.size=7
pdf.y.size=10
rescalingFactor=1.0 #Don't touch this!
printPDF=F
# parFile=readLines(con=parFileName)
########################### Reading par file #########################
parFile=scan(parFileName, character(0), sep = "\n", strip.white = TRUE) # separate each line
pdfFileName=paste(parFileName, ".pdf", sep="")
#--- Clean input file by removing consecutive separators, and keep ---------------
#--- Function to remove separators within a string
removeTrailingSep=function(string, sep) {
temp=strsplit(string, split=sep)
temp2=temp[[1]][nchar(temp[[1]])>0]
cleanStr=temp2[1]
if (length(temp2)>1) {
for (i in 2:length(temp2)) {
cleanStr=paste(cleanStr, temp2[i], sep=sep)
}
}
cleanStr
}
#--- Replace Keep by -9999
replaceKeep=function(string) {
if (grepl("keep", string)) {
return(gsub("keep", "-9999", string))
}
return(string)
}
#Remove both multiple consecutive whitespace and tabs and replace keep keyword
for (i in 1:length(parFile)) {
parFile[i]=removeTrailingSep(parFile[i], sep='\t')
parFile[i]=removeTrailingSep(parFile[i], sep=' ')
parFile[i]=replaceKeep(parFile[i])
}
#-------------------------------------------------------------------------------
#--- Get number of samples on line 2 -----
l.numsamples=parFile[2]
sp.l=unlist(strsplit(l.numsamples, split=' '))
tab.l=unlist(strsplit(l.numsamples, split='\t'))
if (length(sp.l)>=length(tab.l)) {
numSamples=as.numeric(sp.l[1])
separator=" "
} else {
numSamples=as.numeric(tab.l[1])
separator="\t"
}
#--- Reading numbers on separate lines -----
getNumbers=function(start, parFile, numSamples) {
for (i in 1:numSamples) {
curnum=as.numeric(unlist(strsplit(parFile[start+i], split=separator))[1])
if (i==1) {
num=curnum
} else {
num=c(num, curnum)
}
}
num
}
#--- Get population sizes -----------------
start=3
popSizes=getNumbers(start, parFile, numSamples)
#Rescaling pop sizes
popSizes=round(popSizes*rescalingFactor, digits=0)
iniPopSizes=popSizes
#--- Get sample sizes -----------------
readSampleSizesTimesAndInbreedingLevel=function(start, parFile, numSamples) {
for (i in 1:numSamples) {
curLine=unlist(strsplit(parFile[start+i], split=separator))
curSampSize=as.numeric(curLine[1])
curSampTime=0
curInbreeding=0
if (length(curLine)>1) curSampTime=as.numeric(curLine[2])
if (length(curLine)>2) curInbreeding=as.numeric(curLine[3])
if (i==1) {
sampSize=curSampSize
sampTime=curSampTime
inbreeding=curInbreeding
} else {
sampSize=c(sampSize,curSampSize)
sampTime=c(sampTime,curSampTime)
inbreeding=c(inbreeding,curInbreeding)
}
}
list(ss=sampSize, st=sampTime, inb=inbreeding)
}
start=start+numSamples+1
# sampSizes=getNumbers(start, parFile, numSamples)
sampSizesStats=readSampleSizesTimesAndInbreedingLevel(start, parFile, numSamples)
#Rescaling sample times
sampSizesStats$st=round(sampSizesStats$st*rescalingFactor, digits=0)
sampSizes=sampSizesStats$ss
sampTimes=sampSizesStats$st
inbrCoeff=sampSizesStats$inb
#--- Get growth rates -----------------
start=start+numSamples+1
growthRatesInitial=getNumbers(start, parFile, numSamples)
# save this into growthRates which will be used and updated when printing historical events
growthRates=growthRatesInitial
#--- Get number of migration matrices -----------------
start=start+numSamples+1
numMigMat=as.numeric(unlist(strsplit(parFile[start+1], split=separator))[1])
#--- Read migration matrix
readMigMat=function(start, parFile, numSamples) {
for (i in 1:numSamples) {
curmigs=as.numeric(unlist(strsplit(parFile[start+i], split=separator)))
if (i==1) {
migs=curmigs
} else {
migs=rbind(migs, curmigs)
}
}
rownames(migs)=1:numSamples
migs
}
#--- Get migration matrices as a list --------------
start=start+2
migMats=list()
if (numMigMat>0) {
for (i in 1:numMigMat) {
curMigMat=readMigMat(start, parFile, numSamples)
migMats[[i]]=curMigMat
start=start+numSamples+1
}
}
#Rescaling migration rates
if (numMigMat>0) {
for (i in 1:numMigMat) {
migMats[[i]]=migMats[[i]]/rescalingFactor
}
}
#--- Get number of historical events
start=start+1
numHistEvents=as.numeric(unlist(strsplit(parFile[start], split=separator))[1])
###################### HISTORICAL EVENTS HANDLING ##############################
#..... Read Historical Event .......
last.he.time=0
if (numHistEvents>0) {
for (i in 1:numHistEvents) {
start=start+1
#Take care of nomig keyword
nomig=F
if (grepl("nomig", parFile[start])) {
nomig=T
gsub("nomig", "", parFile[start])
}
curHE=as.numeric(unlist(strsplit(parFile[start], split=separator)))
if (nomig) curHE[7]=-1
#Rescaling time of event
curHE[1]=round(curHE[1]*rescalingFactor, digits=0)
if (i==1) {
histEvents=curHE
last.he.time=curHE[1]
} else {
histEvents=rbind(histEvents, curHE)
if (histEvents[i,1] > last.he.time) last.he.time=histEvents[i,1]
}
}
if (numHistEvents>1) rownames(histEvents)=1:numHistEvents
}
names(last.he.time)=""
last.he.time=as.numeric(last.he.time)
yTimeLimit=0
if (last.he.time!=0) yTimeLimit=last.he.time*(1+propLastsegment) #Add propLastsegment to y axis after last event (to draw stuff)
#Reorder events by their times
if (numHistEvents>1) histEvents=histEvents[order(histEvents[,1],decreasing=FALSE),] else {
if (numHistEvents==1) histEvents=matrix(histEvents, nrow=1, byrow = T)
}
endReadParFile=start
############################## PLOTTING THE MODELED SCENARIO ######################################
#--- Graphical functions ...........................................................................
fullHeadArrow=function(x0, y0, x1, y1, length, angle, color="black", weight=1) {
arrows(x0, y0, x1, y1, length, angle, code=2, lty=1, col=color, lwd=weight)
arrows(x0, y0, x1, y1, length, angle*0.80, code=2, lty=1, col=color, lwd=weight)
arrows(x0, y0, x1, y1, length, angle*0.60, code=2, lty=1, col=color, lwd=weight)
arrows(x0, y0, x1, y1, length, angle*0.40, code=2, lty=1, col=color, lwd=weight)
arrows(x0, y0, x1, y1, length, angle*0.20, code=2, lty=1, col=color, lwd=weight)
arrows(x0, y0, x1, y1, length, angle*0.10, code=2, lty=1, col=color, lwd=weight)
}
drawTriangle=function(growth, x, y, size, aspRatio, color) {
if (growth>0) {
x0=x; y0=y
x1=x-size/2; y1=y+size/2*aspRatio
x2=x+size/2; y2=y+size/2*aspRatio
} else {
x0=x-size/2; y0=y
x1=x; y1=y+size/2*aspRatio
x2=x+size/2; y2=y
}
polygon(c(x0, x1, x2, x0), c(y0, y1, y2, y0), col=color)
return(y+size/2*aspRatio)
}
#--- Computing maximum current pop size ............................................................
maxPopSize=popSizes[1]
minPopSize=popSizes[1]
if (length(popSizes)>1){
for (i in 2:length(popSizes)) {
if (popSizes[i]>maxPopSize) maxPopSize=popSizes[i];
if (popSizes[i]<minPopSize) minPopSize=popSizes[i];
}
}
#--- Find min and max pop sizes over the whole population history ..................................
ps=popSizes
isGrowth=FALSE
if (numHistEvents) {
#Need to keep track ofgrowth rates over time
gRates=growthRates
prevTime=0
for (i in 1:numHistEvents) {
he=histEvents[i,]
curTime=he[1]
sink=he[3]+1
resize=he[5]
growth=he[6]
#Begin by resizing pop sizes due to growth since last event
for (j in 1:length(ps)) {
ps[j]=ps[j]*exp(gRates[j]*(curTime-prevTime))
if (ps[j]>maxPopSize) maxPopSize=ps[j]
if (ps[j]<minPopSize) minPopSize=ps[j]
if (ps[j]==0) print(paste("Deme ", j-1, " reaches size zero at time ", curTime, " due to large negative growth (", gRates[j], ")", sep=""))
if (is.infinite(ps[j])) print(paste("Deme ", j-1, " reaches infinite size at time ", curTime, " due to large positive growth (", gRates[j], ")", sep=""))
}
prevTime=curTime
#Handle transformed keep keyword
if(growth!=-9999) {
gRates[sink]=growth
}
#Implement resize
ps[sink]=ps[sink]*resize
if (ps[sink]>maxPopSize) maxPopSize=ps[sink]
if (ps[sink]<minPopSize) minPopSize=ps[sink]
}
}
#-- Function to compute the circle radius for a given pop size ....................................
interpolRadius=function(curSize, minSize, maxSize, minRadius, maxRadius, logScale) {
if(logScale) {
minSize=log10(minSize)
maxSize=log10(maxSize)
curSize=log10(curSize)
}
curRadius=minRadius+(curSize-minSize)*(maxRadius-minRadius)/(maxSize-minSize)
curRadius
}
#============================== ========================
#============================== BEGIN MAIN PLOT ========================
#============================== ========================
library("plotrix")
library("diagram")
if (printPDF) {
pdf(pdfFileName, width=pdf.x.size, height=pdf.y.size)
}
par(xpd=F, mar=c(6,4,3,0.5))
maxRadius=maxRadius*(numSamples+2)
minRadius=minRadius*(numSamples+2)
title=parFileName
plot(x=1:numSamples, type="n", xlab="", ylab="time (gen)", xlim=c(-0.5, numSamples+1.5), cex.main=0.8,
ylim=c(0, yTimeLimit), main=title, xaxt = 'n', log=logScaleAxis, cex.axis=0.9, cex.lab=0.9)
axis(side=1, labels=c("Mig Mat",0:(numSamples-1), " \nTimes"), at=0:(numSamples+1), cex.axis=0.8)
mtext("demes", side=1,line=2, cex=0.8)
w <- par("pin")[1]/diff(par("usr")[1:2])
h <- par("pin")[2]/diff(par("usr")[3:4])
aspRatio <- w/h
# aspRatio=yTimeLimit/(numSamples+2)
#--- Plotting population circles according to their current size and sampling time
for (i in 1:numSamples) {
curRadius=interpolRadius(popSizes[i], minPopSize, maxPopSize, minRadius, maxRadius, drawLogPopSize)
# print(curRadius)
if (i==1) topOfCircle=sampTimes[i]+curRadius*aspRatio else topOfCircle=c(topOfCircle, sampTimes[i]+curRadius*aspRatio)
curColor=popCol
if (sampSizes[i]==0) curColor="white"
if (inbrCoeff[i]==0) {
draw.circle(i, sampTimes[i], radius = curRadius, col=curColor, border=popBorderCol)
} else {
floating.pie(i, sampTimes[i], c(inbrCoeff[i], 1-inbrCoeff[i]), radius = curRadius,
col=c(inbreedColor,curColor), startpos=pi, border=popBorderCol)
}
if (sampTimes[i]>0) {
text(i, sampTimes[i]-curRadius*aspRatio, labels=sampTimes[i], cex=timeProp, col=oldSampColor, pos=1)
}
#Draw a vertical arrow in case of pop growth
curGrowthRate=growthRates[i]
#--- Draw gtrowing or shrinking triangle on top of pop circle to reflect growth type
if (curGrowthRate!=0) {
arLength=0.15
topOfCircle[i]=drawTriangle(curGrowthRate, x=i, y=topOfCircle[i], arLength, aspRatio, color=growthCol)
}
}
#--- Handle first migration matrix .......................
curMigMatNum=0
curvature=0.0075*last.he.time
text(0-migrOffset, 0, labels=0, cex=migMatNameProp, col=migrMatCol)
if (numMigMat) {
curMigMat=migMats[1][[1]]
for (sink in 1:numSamples) {
for (sourc in 1:numSamples) {
if (sink!=sourc & curMigMat[sourc, sink]>0) {
differ=sourc-sink
curvedarrow(from=c(sourc, 0), to=c(sink,0), curve=curvature*(abs(differ)*0.55^abs(differ)), arr.adj=1,
arr.pos=0.5, arr.type="triangle", arr.col=migrMatCol, lwd=1, lty=curvedArrowLTY,
lcol=migrMatCol, arr.length=arrowLength)
if (plotMigrRates) {
curNm=round(curMigMat[sourc, sink]*popSizes[sourc], digits=2)
if (differ>0) {
text(sink+abs(differ)/2, aspRatio*0.15*abs(differ), labels=curNm, cex=migrRateTextSizeCEX, col=migrMatCol)
} else {
text(sourc+abs(differ)/2, -aspRatio*0.15*abs(differ), labels=curNm, cex=migrRateTextSizeCEX, col=migrMatCol)
}
}
}
}
}
}
#---Draw all events on the population tree
lastTime=0
activePops=1:numSamples
numActivePops=numSamples
lastSink=-1
if (numHistEvents) {
for (i in 1:numHistEvents) {
#Extract historical event
he=histEvents[i,]
he.time=he[1]
he.source=he[2]+1 #+1 is due to the use of base 0 for deme number in fsc
he.sink=he[3]+1
he.migr=he[4]
he.resize=he[5]
he.growth=he[6]
if(he.growth==-9999) { #handle transformed keep keyword
he.growth=growthRates[he.sink] #Keep current growth rate
}
he.migrMat=he[7]
if(he.migrMat==-9999) { #handle transformed keep keyword
he.migrMat=curMigMatNum
}
#-- Draw event time
if (i%%2) {
slide=timeOffset
}
else {
slide=-timeOffset
}
#Draw all vertical segments ....................
for (p in 1:length(activePops)) {
#Handle population with older sampling times
if (sampTimes[activePops[p]]>0) {
minTime=max(lastTime, topOfCircle[activePops[p]]);
} else {
minTime=lastTime;
}
if (he.time>sampTimes[activePops[p]] & he.time>topOfCircle[activePops[p]]) {
if (activePops[p]==lastSink | i==1) {
segments(activePops[p], topOfCircle[activePops[p]], activePops[p], he.time)
} else {
segments(activePops[p], minTime, activePops[p], he.time)
}
}
}
#Handle growth rate changes since last event ......................
#Update pop sizes according to current growth rates
if (i>1) {
prev.he=histEvents[i-1,]
branchLength=he.time-lastTime
for (p in 1:length(activePops)) {
curPop=activePops[p]
popSizes[curPop]=popSizes[curPop]*exp(growthRates[curPop]*branchLength);
}
}
#Update growth rate
growthRates[he.sink]=he.growth
lastTime=he.time
#Handle resize of sink pop ........................
popSizes[he.sink]=popSizes[he.sink]*he.resize
curRadius=interpolRadius(popSizes[he.sink], minPopSize, maxPopSize, minRadius, maxRadius, drawLogPopSize)
# popRadius[he.sink]=curRadius
topOfCircle[he.sink]=he.time+curRadius*aspRatio
draw.circle(he.sink, he.time, radius = curRadius, col=popCol, border=popBorderCol)
#--- Draw growing or shrinking triangle on top of pop circle to reflect growth type
if (he.growth!=0) {
arLength=0.15
topOfCircle[he.sink]=drawTriangle(he.growth, x=he.sink, y=topOfCircle[he.sink], arLength, aspRatio, color=growthCol)
}
#Handle population fusion .........................
if (he.migr>=1 & he.sink!=he.source) { #This is a population fusion
if (numActivePops==numSamples) removedPops=he.source else removedPops=c(removedPops, he.source)
numActivePops=numActivePops-1
activePops=(1:numSamples)[-removedPops]
#Draw connecting arrows from source to sink
fullHeadArrow(he.source, he.time, he.sink, he.time, length=0.15, angle=20)
#Redraw time with the right color
text(numSamples+1+slide, he.time, labels=he.time, cex=timeProp, col=popFusionColor)
} else {
#--- Handle admixture event ........................
if (he.migr>0 & he.migr<1) {
#Draw connecting arrows from source to sink
segments(he.source, he.time, he.sink, he.time, col=admixCol, lty=2)
if (he.sink>he.source) {
fullHeadArrow(he.sink-0.15, he.time, he.sink, he.time, length=0.15, angle=20, color=admixCol)
} else {
fullHeadArrow(he.sink+0.15, he.time, he.sink, he.time, length=0.15, angle=20, color=admixCol)
}
#Redraw time with the right color
text(numSamples+1+slide, he.time, labels=he.time, cex=timeProp, col=admixCol)
}
else text(numSamples+1+slide, he.time, labels=he.time, cex=timeProp, col=timeCol)
}
#--- Handle migmat change ............................
if (i!=numHistEvents) nextTime=histEvents[i+1,1] else {
nextTime=yTimeLimit
}
time2DrawArrows=(he.time+nextTime)/2
if (he.migrMat!=curMigMatNum) {
if (he.migrMat>-1) {
curMigMat=migMats[he.migrMat+1][[1]]
for (sink in 1:numSamples) {
for (sourc in 1:numSamples) {
if (sink!=sourc & curMigMat[sourc, sink]>0) {
differ=sourc-sink
curvedarrow(from=c(sourc, time2DrawArrows), to=c(sink, time2DrawArrows), curve=curvature*(abs(differ)*0.55^abs(differ)),
arr.adj=1, arr.pos=0.5, arr.type="triangle", arr.col=migrMatCol, lwd=1, lty=curvedArrowLTY,
lcol=migrMatCol, arr.length=arrowLength)
if (plotMigrRates) {
#Write Nm values
curNm=round(curMigMat[sourc, sink]*popSizes[sourc], digits=2)
if (differ>0) {
text(sink+abs(differ)/2, time2DrawArrows+aspRatio*0.15*abs(differ), labels=curNm, cex=migrRateTextSizeCEX, col=migrMatCol)
} else {
text(sourc+abs(differ)/2, time2DrawArrows-aspRatio*0.15*abs(differ), labels=curNm, cex=migrRateTextSizeCEX, col=migrMatCol)
}
}
}
}
# curMigMatNum=he.migrMat
}
}
#--- Draw separation between migration matrices numbers on the left
segments(-migMatLineLength/2, he.time, migMatLineLength/2, he.time, lty=3, col=migrMatCol)
}
curMigMatNum=he.migrMat
#Output current valid migration matrix
if (i%%2) {
slide=migrOffset
}
else {
slide=-migrOffset
}
if (he.migrMat>-1) {
migText=he.migrMat
curCex=migMatNameProp
} else {
migText="nomig"
curCex=migMatNameProp/2
slide=slide*2
}
text(0+slide, time2DrawArrows, labels=migText, cex=curCex, col=migrMatCol)
lastSink=he.sink
}
}
#-- Draw last branch
segments(activePops[1],topOfCircle[activePops[1]], activePops[1], yTimeLimit)
#============================== PLOT LEGENDS IN MARGINS ========================
#Compute space available in margin
minY.coo=grconvertY(0, from="nic", to="user")
par(xpd=NA)
#--- Draw population size scale with circles of different sizes .................
maxOrder=ceiling(log10(maxPopSize))
minOrder=floor(log10(minPopSize))
popSizeRadius=10^(maxOrder:minOrder)
winWidth=numSamples+2
ypos=3/4*minY.coo
text(x=-winWidth/10*1.2, y=ypos, labels="Pop. \nsizes ", cex=.8, pos=2)
for (i in 1:length(popSizeRadius)) {
curRadius=interpolRadius(popSizeRadius[i], minPopSize, maxPopSize, minRadius, maxRadius, drawLogPopSize)
# print(curRadius)
if (curRadius>0) {
xpos=-winWidth/10+(i-1)*winWidth/10
draw.circle(xpos, ypos, radius=curRadius, col=popCol, border=popBorderCol)
}
text(xpos, ypos-abs(ypos)*0.1, popSizeRadius[i], cex=0.7, pos=1, col="black")
}
#--- Legend for growing or shrinking populations ...............................
if (isGrowth) {
x=winWidth-1.5*winWidth/10; y=ypos+abs(ypos)*0.2
text(x, y-abs(ypos)*0.1, labels="Populations", cex=0.8)
x=winWidth-2*winWidth/10; y=ypos-abs(ypos)*0.1
drawTriangle(1, x, y, size=0.15, aspRatio, color=growthCol)
text(x, y-abs(ypos)*0.1, labels="growing", pos=NULL, cex=0.7)
x=winWidth-winWidth/10
drawTriangle(-1, x, y, size=0.15, aspRatio, color=growthCol)
text(x, y-abs(ypos)*0.1, labels="shrinking", pos=NULL, cex=0.7)
}
#if (printPDF) dev.off()
dev.off()
## null device
## 1